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ABSTRACT 

Growth of massive black holes (MBHs) in galactic centers comes mainly from gas accretion dur- 
ing their QSO/AGN phases. In this paper we apply an extended Soltan argument, connecting the 
local MBH mass function with the time-integral of the QSO luminosity function, to the demogra- 
phy of MBHs and QSOs from recent optical and X-ray surveys, and obtain robust constraints on 
the luminosity evolution (or mass growth history) of individual QSOs (or MBHs). We find that the 
luminosity evolution probably involves two phases: an initial exponentially increasing phase set by 
the Eddington limit and a following phase in which the luminosity declines with time as a power law 

(with a slope of ~ —1.2 1.3) set by a self-similar long-term evolution of disk accretion. Neither 

an evolution involving only the increasing phase with a single Eddington ratio nor an exponentially 
declining pattern in the second phase is likely. The period of a QSO radiating at a luminosity higher 
than 10% of its peak value is about 2-3xl0 8 yr, during which the MBH obtains ~ 80% of its mass. 
The mass-to-energy conversion efficiency is ~ 0.16±0.04lg' 05 , with the latter error accounting for the 
maximum uncertainty due to Compton-thick AGNs. The expected Eddington ratios in QSOs from 
the constrained luminosity evolution cluster around a single value close to 0.5-1 for high-luminosity 
QSOs and extend to a wide range of lower values for low-luminosity ones. The Eddington ratios for 
high luminosity QSOs appear to conflict with those estimated from observations (~ 0.25) by using 
some virial mass estimators for MBHs in QSOs unless the estimators systematically over-estimate 
MBH masses by a factor of 2-4. We also infer the fraction of optically obscured QSOs ~ 60 — 80%. 
The constraints obtained above are not affected significantly by MBH mergers and multiple-times of 
nuclear activity (e.g., triggered by multiple times of galaxy wet major mergers) in the MBH growth 
history. We discuss further applications of the luminosity evolution of individual QSOs to obtain- 
ing the MBH mass function at high redshifts and the cosmic evolution of triggering rates of nuclear 
activity. 

Subject headings: black hole physics - galaxies: active - galaxies: evolution - galaxies: nuclei - quasars: 
general - cosmology: miscellaneous 



1. INTRODUCTION 

Massive black holes (MBHs), probably remnants of 
QSOs (Lynden-Bell 1969), have been detected in the nu- 
clei of many nearby galaxies (Kormendy & Richstone 
1995; Magorrian et al. 1998; Richstone et al. 1998; Kor- 
mendy & Gebhardt 2001; Ferrarese & Ford 2005). How 
do these local MBHs form and evolve, and what is the 
most important mechanism shaping the mass distribu- 
tion of MBHs? The current consensus is that the lo- 
cal MBHs obtained their mass mainly through accretion 
during phases of nuclear activity when they appeared as 
QSOs/ AGNs, 4 similar to the ones seen now in the dis- 
tant universe (e.g., Yu & Tremaine 2002; Yu & Lu 2004a; 
Marconi et al. 2004; Shankar et al. 2004; Barger et al. 
2005; Hopkins et al. 2006; Shankar et al. 2007). The evo- 
lution of mass accretion onto a MBH is equivalent to the 
luminosity evolution, given the mass-to-energy conver- 
sion efficiency, and is recorded in the luminosity function 
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(LF) of QSOs. However, the QSO LF depends mainly 
on two functions: (1) Q(z; M,^), the rate of nuclear ac- 
tivity triggered at different redshifts z for MBHs with 
present-day mass M. i0 ; (2) £(r;M #! o), the luminosity 
evolution history of a QSO, of which the remnant MBH 
has a present-day mass M.,q, as a function of the age of 
its nuclear activity r. One cannot derive these two func- 
tions only from the knowledge of the QSO LF without 
additional assumptions. 

In an extended version of the Soltan (1982) argument, 
the local MBH mass distribution function (BHMF) is 
related to QSOs found in the distant universe by the 
simple integral equation 



riM.{M,fi, to) x 



T lt (M. :0 )P(L\M. :0 )dM..o, (1) 

where to is the present cosmic time, um. ^o) is the 

local BHMF, defined so that um. (-^.,0, t )dM,^ gives 
the number density of local MBHs with present-day mass 
in the range M., — > M. j0 + dM, y o, ^l(L, z) is the QSO 
LF, defined so that \& l(L, z)dL gives the comoving num- 
ber density of QSOs with nuclear luminosity in the range 
L — > L + dL at redshift z, 
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is the time interval (or the QSO lifetime) in which that 
a MBH with present-day mass M #i o appeared as a QSO, 
and Tk(L,M, t o) (k = 1,2, ...) are the roots of the equa- 
tion £(t; M. j0 ) — L = (see details of the derivation in Yu 
& Lu 2004a). Here £(r; M. i0 ) represents the luminosity 
of a QSO and its associated MBH with present-day mass 
M. o at a time t after the triggering of nuclear activ- 
ity. The value of Tit depends on the detailed definition of 
"active nuclei" or the lower threshold set to the nuclear 
luminosity. Finally, 
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is the probability distribution function of the nuclear 
(bolometric) luminosity L over the growth history of the 
MBH. The right-hand-side of equation (1) gives the total 
time spent per unit L at luminosity L by the progeni- 
tors of all the local MBHs in a unit comoving volume, 
which should be the time integral of the QSO LF, i.e., 
the left-hand-side of the equation. Multiplying equation 
(1) by the BH mass accretion rate (1 — e)L/(ec 2 ) = M, 
(see eqs. 26 and 27 below), where e is the mass-to-energy 
conversion efficiency and c is the speed of light, and then 
integrating it over cosmic time t reduces to the Soltan 
(1982) argument (Yu & Lu 2004a). Provided that two 
basic quantities, i.e., the local BHMF and the QSO LF, 
can be observationally determined with sufficient accu- 
racy, the kernel T\t(M,fi)P(L\M,^), containing informa- 
tion on the luminosity evolution history of individual 
QSOs/MBHs, may be solved from the integral equa- 
tion (1). Therefore, the extended Soltan argument is 
expected to give robust but more detailed constraints on 
the growth of MBHs than the simple energetic argument 
due to Soltan (1982). 

As an alternative approach to the theoretical mod- 
els based on the hierarchical co-evolution of MBHs and 
galaxies/galactic halos studied intensively in the liter- 
ature (e.g., Efstathiou & Rees 1988; Haehnelt & Rees 
1993; Haehnelt et al. 1998; Kauffmann & Haehnelt 2000; 
Granato et al. 2001; Wyithe & Loeb 2003; Volonteri et 
al. 2003; Croton et al. 2006; Bower et al. 2006; Malbon 
et al. 2007), in this paper we use the integral equation 
(1) to statistically constrain the growth history of indi- 
vidual MBHs or C. The advantages of this approach are: 
(1) the accretion history of individual QSOs, £(t; M. o), 
is isolated from the triggering rate of nuclear activity, 
Q(z;M,fi), which is presumably associated with merg- 
ers of galaxies or instabilities of galactic disks; and (2) 
it is free of the many adjustable parameters introduced 
in the co-evolution models and probably also avoids un- 
certain assumptions on seed BHs. Note that these two 
functions, Q(z\M t ^) and £(t;M. i0 ), are mixed in the 
differential continuity equation for BHMF evolution pre- 
sented in Small & Blandford (1992, see also Cavaliere 
et al. 1971, Cavaliere & Padovani 1989, and Caditz & 
Petrosian 1990), which is widely used in studying the 
growth of MBHs (e.g., Marconi et al. 2004; Shankar et 
al. 2007). Using the luminosity evolution curves, i.e., 
£(t;M. j0 ), obtained from numerical simulations of col- 
liding galaxies, Hopkins et al. (2006) elaborated a unified 
model for the origin of QSOs and MBHs (see also their 
other papers listed therein) . A possible concern with that 



approach is that simulations of colliding galaxies have a 
spatial resolution much larger than the scale of accretion 
disks around MBHs and therefore may not reflect the 
real luminosity evolution, as the disk accretion is prob- 
ably self-regulated in the vicinity of MBHs rather than 
being directly determined by the material infall rate from 
a much larger scale or the Bondi-accretion rate (see dis- 
cussions in § 4). (For another model of the possible light 
curve, see Ciotti & Ostriker 2007.) 

Estimating the local BHMF can be done with recent 
advances in observations (e.g., Salucci et al. 1999; Aller 
& Richstone 2002; Yu & Lu 2004a; Marconi et al. 2004; 
Shankar et al. 2004; Lauer et al. 2007a; Tundo et al. 
2007). First, MBHs arc believed to exist in the nuclei of 
most, if not all, nearby galaxies (Kormendy & Richstone 
1995; Magorrian et al. 1998; Richstone et al. 1998; Kor- 
mendy & Gebhardt 2001; Ferrarese & Ford 2005). Sec- 
ond, it has been well established that tight correlations 
exist between the MBH mass and various galactic proper- 
ties, such as mass, luminosity, stellar velocity-dispersion, 
light concentration and binding energy of the hot com- 
ponents of galaxies (here hot components mean either 
ellipticals or spiral bulges; Kormendy & Richstone 1995; 
Magorrian et al. 1998; Ferrarese & Merritt 2000; Geb- 
hardt et al. 2000; Tremaine et al. 2002; Haring & Rix 
2004; Marconi & Hunt 2004; Graham et al. 2001; Aller 
& Richstone 2007). Third, the luminosity or velocity- 
dispersion functions of nearby galaxies have been well 
determined by large surveys such as the Sloan Digital Sky 
Survey (SDSS; Blanton et al. 2003; Bernardi et al. 2003; 
Shcth et al. 2003). Combining the correlation between 
the MBH mass and galaxy velocity dispersion (or lumi- 
nosity) with the velocity-dispersion (or luminosity) dis- 
tribution of nearby galaxies, we estimate the local BHMF 
in§ 2. 

In the past several years, the QSO LF has been deter- 
mined over unprecedentedly large luminosity and redshift 
ranges both from optical surveys such as the Two Degree 
Field QSO Redshift Survey (2Qz) and SDSS, and from 
X-ray surveys by ASCA, Chandra and XMM-Newton. 
For example, the optical QSO LF has been obtained 
over the redshift range 0.4 < z < 2.1 and the magni- 
tude range Mb, < —22.5 using a sample of more than 
15,000 QSOs from 2Qz (Croom et al. 2004); Richards 
et al. (2006a) estimated the QSO LF over a larger red- 
shift range (0.3 < z < 5), but only for bright QSOs, 
using a homogeneous statistical sample of 15,343 QSOs 
drawn from SDSS Data Release 3; using the COMBO-17 
data, Wolf et al. (2003) estimated the LF for faint QSOs 
over the range 1.2 < z < 4.8; and Jiang et al. (2006) 
estimated the QSO LF over the range 0.5 < z < 3.6 by 
using a deep survey of faint QSOs in the SDSS. Obscured 
(or type 2) QSOs may be missed in the optical surveys 
but can be detected in hard X-ray surveys. La Franca et 
al. (2005) use 508 AGNs to estimate the hard X-ray LF 
(HXLF; 2-10 keV) over the range < z < 2.5 by com- 
bining data from XMM-Newton (Lockman hole) and the 
Chandra Deep Field (CDF). Barger et al. (2005) use a 
spectroscopically complete deep and wide-area Chandra 
survey to estimate the HXLF (2 — 8 keV) over the range 
< z < 5. Silverman et al. (2008) measure the HXLF 
(2 — 8 keV) up to z ~ 5 with fewer uncertainties by com- 
bining the observations from the CDF and the Chandra 
Multiwavclcngth Project. Combining all these observa- 
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tions, the time integrals of the QSO LF are estimated in 
§3. 

In § 4, we assume several models for the luminosity 
evolution history of individual QSOs, i.e., £(t;M. ;0 ), 
and then apply the models and the observational BHMF 
and QSO LF to equation (1) to give constraints on the 
growth of individual MBHs and the associated parame- 
ters, specifically, the efficiency (mainly determined by the 
spin of a MBH), the lifetime of nuclear activity, and the 
long-term evolution of disk accretion etc. We find that 
a reference model for the luminosity evolution history 
of individual QSOs, i.e., an initial rapid accretion phase 
with a rate close to the Eddington limit and then a fol- 
lowing power-law declining phase set by the self-similar 
long-term evolution of disk accretion (M. oc r -7 , and 
7 ~ 1.2 — 1.3), can satisfy the extended Soltan argument 
(eq. 1) well. Using the reference model for £(M,.o, r), we 
discuss the role of obscuration in the BH growth history 
in § 5 and find that obscuration is unlikely to be solely an 
evolutionary effect. The luminosity (or accretion-rate) 
evolution constrained by the extended Soltan argument 
also implies a distribution of Eddington ratios (i.e., the 
accretion rate in units of the Eddington limit) in QSOs. 
In § 6, we particularly discuss its distribution expected 
from the models and compare them with observations. 
In § 7, by using toy models, we discuss the effects of 
BH mergers on our results, which are shown to be in- 
significant. In § 8, we discuss further implications of the 
luminosity evolution obtained from the extended Soltan 
argument. Together with the QSO LF, £(r; M >i0 ) can 
be used to further derive the BHMF at redshift z and 
the triggering rate of nuclear activity Q(z; M, t o). Given 
£(t; M #; o) and Q{z; M. j0 ), many statistical properties of 
QSOs can be inferred and comparison of them with ob- 
servations may further deepen our understanding of the 
growth of MBHs. Conclusions are given in § 9. 

In this paper we set the Hubble constant as Hq = 
100/ikm s _1 Mpc~ ; and if not otherwise specified, the 
cosmological model used is (fl m ,flA,h) — (0.3,0.7,0.7). 

2. THE MASS FUNCTION OF MBHS AT z = 

Studies of central MBHs in nearby galaxies have re- 
vealed strong correlations between the BH mass and the 
velocity dispersion (or luminosity, or other properties) of 
the hot stellar component of the host galaxy (e.g., Ko- 
rmendy & Richstone 1995; Magorrian et al. 1998; Fer- 
rarese & Merritt 2000; Gebhardt et al. 2000; Tremaine 
et al. 2002; Haring & Rix 2004; Marconi & Hunt 2004; 
Graham et al. 2001; Aller & Richstone 2007; Hopkins ct 
al. 2007b). We first present several latest fits of these cor- 
relations (i.e., the M,fi — a relation and the M #j o — -^buige 
relation) and then present the observational velocity- 
dispersion (and luminosity) distribution of nearby galax- 
ies. By combining them, we estimate the local BHMF in 
§ 2.3. 

2.1. The M, t o — a and M. o — -^buigo relations 

Lauer et al. (2007a) show that the logarithm of the BH 
mass at a given velocity dispersion a has a mean value 
given by 

(log M. )0 | logcr) = (8.29 ± 0.07) + (4.13 ± 0.32) x 

1O ^200kn^)' (4) 



which is fitted in the (log M #i o, logcr) space. The mean 
value at a given F-band absolute magnitude My is given 
in the same paper as 

(logM.olMv) = (8.67 ± 0.09) - ( L32 ± °- 14 ) r My + 2 2). 

2.5 

(5) 

The intrinsic scatters around the relations above are not 
reported in Lauer et al. (2007a). (Hereafter the in- 
trinsic scatters in logM.^o are noted as Am. >0 -<t and 
A M.,o-Lbui g c for tnc M *,o - cr relation and M. j0 - L hulgc 
relation, respectively.) Based on the same sample, 
Tremaine et al. (2002) estimate that the intrinsic scatter 
in logM #i o for the M.^ — cr relation, i.e., Am. - CT5 should 
be not larger than 0.25 — 0.3 dex. The latest fit of the 
M. — cr relation by Hu (2008), which is consistent with 
that given by Lauer et al. (2007a) on the zero point and 
the slope, also gives an upper limit to the intrinsic scat- 
ter ~ 0.25 dex. Note also the zero point in equation (4) 
is larger than that obtained by Tremaine et al. (2002) by 
0.10 dex, but roughly consistent with statistical errors. 

The estimates of the M,, — & and M. )0 — ibuige rela- 
tions in Bernardi et al. (2007) are given by: 

(log M. )0 | logcr) = (8.21 ± 0.05) + (3.83 ± 0.10) x 

l0g ( 200km s-0 ' (6) 

and 

(logM., |M r ) = (8.57 ± 0.10) - t 1 - 30 ^ - 10 ) {Mr + 22 ), 

(7) 

with intrinsic scatters not larger than 0.22 ± 0.05 dex 
and 0.33 ± 0.08 dex, respectively. Another set of fits to 
equation (7) by the same authors (Tundo et al. 2007) 
finds a slope 1.30 ±0.15 and the zero point is 8.68 ±0.10, 
consistent with statistical errors. If we convert M r to My 
with M r — My — 0.37 adopted for early-type galaxies 
(Fukugita et al. 1996), we find the zero point in equation 
(7) is larger than that in equation (5) by 0.09 dex. 

The typical difference in the zero point among different 
sets of fits to the M,, — a ( or M,_ — ibuige) relation is 
< 0.10 dex, which is roughly consistent with the statisti- 
cal errors in the zero point estimation. The difference in 
the slope among different sets of fits to the M,^ — a rela- 
tion is quite large compared to the statistical errors in the 
fits, for example, it is 4.13 ±0.32 in Lauer et al. (2007a), 
3.83 ± 0.10 in Bernardi et al. (2007), and 4.86 ± 0.43 in 
Ferrarese & Ford (2005) (for details of the difference in 
the slope see discussions in Tremaine et al. 2002). Note 
that Aller & Richstone (2007) investigate an alternative 
relation to the M #i o — cr relation; and they find that the 
relation between the MBH mass and the bulge gravita- 
tional binding energy is as good as the M..o— cr relation in 
predicting MBH mass but with a slope much more stable 
regarding of changes in the fitting algorithm. A detailed 
study by Novak et al. (2006) demonstrates that the up- 
per limit to the intrinsic scatter is ~ 0.2 — 0.3 dex in the 
M #j o— cr relation and is ~ 0.3—0.4 dex in the M,^— Lbuige 
relation for currently available samples. Below we adopt 
Am,, -» ~ 0.3 dex relation and A M . i0 -z, bulge ~ 0.4 dex 
if not otherwise specified. 

Among the subtle differences in zero points, slopes and 
intrinsic scatters of those relations estimated by different 
groups, the intrinsic scatter would be the most significant 
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one for the purpose of studying MBH growth, because it 
may affect the estimates of the abundance of MBHs at 
the high-mass end (> 10 9 M Q ) by orders of magnitude 
(as shown in Fig. 1 below; see also discussions in Yu 
& Lu 2004a; Marconi et al. 2004; Lauer et al. 2007a; 
Tundo et al. 2007), and this abundance is crucial for our 
understanding of the growth of the most massive BHs in 
bright QSOs. 



2.2. The velocity- dispersion distribution function 
the luminosity function of nearby galaxies 

We define n a (a, t) as the comoving velocity-dispersion 
function of the hot stellar components of local galax- 
ies so that n CT (er, t )dcr represents the number density of 
local galaxies in the range a — > a + da. The velocity- 
dispersion distribution n a {o~, to) includes the contribu- 
tion from both early- type galaxies n%(a,to) and bulges 
of late- type galaxies rv a (cr, to), that is, 



n a (a,to) = <(cr,t ) +n\{o-,to). 



(8) 



• The velocity-dispersion distribution in early-type 
galaxies has been estimated by recent studies of a 
sample of early-type galaxies at z < 0.3 obtained 
by the SDSS (see eq. 4 in Sheth et al. 2003, and 
Bernardi et al. 2003): 



t a ) 



o-y exp[-(o- /**)&] (3 
a J T(a/p) a 1 



(9) 



where the best-fit values of (0», cr*, a, (3) are 
(0.0020+0.0001, 88.8±17.7, 6.5±1.0, 1.93±0.22), 0* 
is the comoving number density of local early-type 
galaxies in units of /iq 7 Mpc -3 , and cr* is in units 
of km s _1 . The brightest cluster galaxies (BCGs) 
are probably under-represented in the above sam- 
ple (Lauer et al. 2007a) . We correct this by adding 
the number density of BCGs to equation (9), where 
the number density of BCGs with a > 350km s _1 is 
estimated from the sample of Bernardi et al. (2006) 
as done in Lauer et al. (2007a). If the scatter in 
the M., — o~ (or M., — Lbuigo) relation is not sig- 
nificantly smaller than 0.3 dex (or 0.4 dcx), this 
correction is not significant because most of high- 
mass MBHs (larger than a few 1O 9 M ) actually 
come from 'modest' galaxies with unusually large 
MBHs for their velocity dispersions or luminosities 
(see the dependence of the BHMF on different val- 
ues of the scatter in Fig. 1; see also Lauer et al. 
2007b). 

• The velocity-dispersion distribution in late-type 
galaxies rt} a (o~, to) may be estimated in the following 
ways, (i) Following Sheth et al. (2003), the LF of 
the late-type galaxies can be obtained by subtract- 
ing the LF of the early-type galaxies (Bernardi et 
al. 2003) from the total LF of all galaxies (Blanton 
et al. 2003). (ii) Following Sheth et al. (2003), the 
distribution of the circular velocity w c in late-type 
galaxies may be obtained by using the LF of the 
late-type galaxies obtained above and the follow- 
ing Tully-Fisher relation (Giovanelli et al. 1997) 



where Mj is the absolute magnitude of the galax- 
ies in the / band, with accounting for the intrinsic 
scatter around relation (10) and the inclination ef- 
fects of galaxies (see details in Sheth et al. 2003). 
(iii) The velocity-dispersion function of late-type 
galaxies can be obtained by using the circular- 
velocity distribution of the late-type galaxies ob- 
tained above and the following relation between the 
circular velocity and the velocity dispersion of the 
bulge component (see eq. 3 in Baes et al. 2003, and 
also Ferrarese 2002): 



log( 



200km s" 



log 



2v r , 



km s _1 



1.00- (M/ -51og/i . 7 )/7.95, (10) 



= (0.96 ±0.11) log. 
v ' & V 200^8-!. 

+ (0.21 ±0.023). (11) 

The intrinsic scatter of relation (11) is small (< 
0.15 dex, see Fig. 1 in Baes et al. 2003) and will be 
ignored in our calculations. We could also simply 
use a = v c /\/3 (e.g., see problem 4.35 in Binney & 
Trcmaine 2008) to estimate cr, which only induces 
a slight difference in estimating the BHMF. Rela- 
tion (11) may not hold for a < 80km s , which 
corresponds to M., < 4 x 10 6 M Q according to the 
M. ;0 — cr relation above (eqs. 4 and 6), but this is be- 
yond the main range which we focus on in § 4. Note 
that the local BHMF for mass M., > 4 x 1O 7 M 
is dominated by the early-type galaxies (see also 
Fig. 1 in Yu & Lu 2004a). 

The LF of galaxies is conventionally described by the 
Schechter (1976) function: 

$(M) = 0.41n(10)(M0-°- 4(M - M * )(Q+1) x 

exphlO- ' 4 ^-^], (12) 

where $(M)cLM gives the comoving number density of 
galaxies with absolute magnitude in the range M — > 
M+dM. Based on observations by the SDSS (Blanton et 
al. 2003), the best fit parameters [0*/(lO~ 2 ^ 7 Mpc~ 3 ), 
M* — 51og/i 7, a] of the LFs are (6.36 ±0.23,' -18.62 ± 
0.02,-0.89±0.03) in the g band and (4.34±0.12,-19.67± 
0.01, — 1.05 ± 0.01) in the r band, respectively. Here M 
is the absolute magnitude of a galaxy (not just of its 
hot stellar component). We can crudely estimate the lu- 
minosity of the hot stellar component of a galaxy, for 
which the relations in equations (5) and (7) are applied, 
from the total luminosity of the galaxy L ga i by setting 
-kbuige = (i g ai/i*)/(l + Lgpi/LJLggx (e.g., Tundo et al. 
2007). With this modification, the BHMF can be esti- 
mated using either the M. ; o — My relation (eq. 5) or the 
M.,o — M r relation (eq. 7) and the galaxy LF in the g 
band (with a color correction of g = My + 0.41; Fukugita 
et al. 1995) or the r band. 

2.3. n M . (M.,o, t ) 
We show in Figure 1 the BHMF obtained by combin- 
ing the M.,o — cr (or M. o — ^tmige) relation with the 
velocity-dispersion (or luminosity) distribution function 
of local galaxies (e.g., see eq. 44 in Yu & Lu 2004a). 
Our calculations show that the uncertainties in the in- 
trinsic scatter of the M., — o (or Af.,o — ^buigo) relation 
may affect estimates of the BHMF significantly at the 
high-mass end (see Fig. 1). To illustrate this effect, we 
assume that the intrinsic scatters in the M., — cr (or 
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M.,o — ibuige) relation by Lauer et al. (2007a) (eqs. 4 
and 5) and by Bernardi et al. (2007) (eqs. 6 and 7) are 0, 
0.2 and 0.3 dex (or 0, 0.3 and 0.4 dex), respectively. With 
the intrinsic scatter of the M #; o — cr (or M.,o — Lbuigc) re- 
lation ~ 0.3 dex (or <~ 0.4 dex), the estimated abundance 
of MBHs at the high-mass end (> 10 9 M Q ) is larger than 
that estimated from a zero intrinsic scatter by orders of 
magnitude (see the upper panels of Fig. 1). The differ- 
ence in the slope and the zero point among different sets 
of fits to the M. ; o — cr (or M.,o — Lbulge) relation may 
also affect the estimates of the abundance of MBHs at 
the high-mass end, but its effects arc substantially less 
significant compared to that of the intrinsic scatter (see 
Fig. 1 and also Yu & Lu 2004a). 

As shown in the bottom left panel of Figure 1, the 
abundance of MBHs estimated from the M., — ibuige 
relation is larger than that from the M.,o — cr relation 
roughly by a factor ~ 2 if both relations are adopted 
from Lauer et al. (2007a) (see also discussions in Lauer 
et al. 2007a and Tundo et al. 2007), but the shapes are 
similar. One possible reason for this discrepancy in abun- 
dance is that the local MBH sample used to derive the 
M.. 

o — ibuigo relation is biased relative to the SDSS 
galaxy sample as discussed in Yu & Tremaine (2002) 
and Bernardi et al. (2007). (The other possibility is sys- 
tematic differences in measurements of luminosity or ve- 
locity dispersion between other surveys and the SDSS.) 
If we correct this 'bias' with the recipe introduced in 
Tundo et al. (2007), then the BHMF estimated from the 
M.. o — ibuigo relation is almost the same as that esti- 
mated from the M #; o — cr relation at the high-mass end 
(M..o si a few 10 8 M Q ), as shown in the bottom right 
panel of Figure 1 . The remaining discrepancy at the low- 
mass end is possibly due to uncertainties in the estima- 
tion of the bulge luminosity from the total luminosity for 
late-type galaxies. For example, recent studies by Lau- 
rikaincn et al. (2005) and Graham & Worlcy (2008) have 
shown that the bulge-to-total luminosity ratio (B/T ra- 
tio) is around 0.24 for SO galaxies, which is substantially 
smaller than the previous estimates (~ 0.6; e.g., Fukugita 
et al. 1998). According to these new estimates, the B/T 
ratio adopted in Tundo et al. (2007) may be an overesti- 
mate at least for SO galaxies, and thus the BHMF at the 
low-mass end < 10 8 Mq is probably substantially overes- 
timated. (The B/T ratio adopted in other estimates of 
the BHMF may be also overestimated; e.g., Marconi et 
al. 2004.) It is anticipated that the BHMF at the low- 
mass end estimated by using the M..o — itmige relation 
will be closer to that estimated by using the M., — cr 
relation if adopting a more realistic B/T ratio for spiral 
galaxies. In § 4, we adopt the BHMF obtained from the 
M. o — a relation given by Lauer et al. (2007a) with an 
intrinsic scatter of 0.3 dex as the reference BHMF, if not 
otherwise specified. 

In addition to the uncertainty on the local BHMF due 
to the intrinsic scatter in the M.,o — cr (or M. ; o — Lbuigc) 
relation, the local BHMF suffers other uncertainties, in 
particular, the uncertainties in estimating the M. o — & 
(or M., — ibuigo) relation [e.g., due to (1) limited mass 
range and small samples; (2) being restricted to ellipti- 
cals, and little is known about late-type galaxies; (3) de- 
termining M.^ is difficult and may be underestimated, 
especially for BCGs] and the uncertainties in estimating 



TABLE 1 

The total mass density of massive black holes 



Method 


Reference 


Note 


P«,0 

10 s M 


M.,o - a 


Lauer07a 




q O + 0.7 

3-O-0.6 


M.,o - a 


Bernardi07 




O q + 0.5 

O.O_ .4 


M.,o - a 


FF05 




o fi +0.7 

^•"-o.e 


M.,o — Lbuigc 


Lauer07a 




7.61?;? 


M.,0 — Lbulge 


Bernardi07 






M.,0 — Lbulge 


Lauer07a 


Bias corrected 


4-9±i; 


M.,0 — Lbulge 


Bernardi07 


Bias corrected 





Note. — The total mass density of massive black holes esti- 
mated from the M.,o — cr (or M.,o — Lbulge) relation obtained by 
different authors. The references for the M,fi—a (or M.,o — Lbui gc ) 
relation are listed in Column 2, and Laucr07a, Bcrnardi07, FF05 
represent Lauer et al. (2007a), Bernardi et al. (2007) and Ferrarese 
& Ford (2005), respectively. 

the velocity-dispersion (or bulge luminosity) distribution 
in late-type galaxies. 

The total mass density of local MBHs can be estimated 
from the BHMF. The differences in the zero point, the 
slope and the intrinsic scatter among the relations esti- 
mated by different groups could cause at most a 20-30% 
difference in the total mass density of local MBHs (as 
shown in § 2.1). For example, adopting the M., — cr 
relation given by Lauer et al. (2007a) yields a total mass 
density of MBHs ~ 3.81°;^ x lO 5 /^ . 7 M Mpc~ 3 , which is 
larger than that obtained by Yu & Tremaine (2002) by a 
factor of <~ 1.3 mainly due to the larger zero point of the 
M.,o — cr relation in Lauer et al. (2007a) adopted here. 
We show in Table 1 a few estimates of the total mass 
density of local MBHs obtained from the M. o — cr (or 
M. o — Lbulge) relation given by different authors. The 
errors are obtained by accounting for the uncertainties in 
the M.,o — cr (or M., — Lbuigc) relation and the galaxy 
velocity-dispersion (or the luminosity) distribution func- 
tion. (For other estimates of the total mass density of 
local MBHs, see Tab. 3 in Graham 2007.) The total 
mass density obtained from the M.,o — Lbulge relation 
is about a factor of ~ 2 larger than that obtained from 
the M.,o — cr relation, which is consistent with that in 
Yu & Tremaine (2002) (see also discussions for the rea- 
sons of this discrepancy in Tundo et al. 2007). If we 
use the recipe introduced by Tundo et al. (2007) to cor- 
rect the possible bias in MBH masses estimated from the 
M. o — -lbuigc relation, the corrected total mass densities 
are still larger than that obtained from the M. o — cr rela- 
tion but now appears to be consistent within statistical 
errors (see Tab. 1). Furthermore, considering that the 
B/T ratio for spiral galaxies adopted in the estimates of 
total BH mass density using the M.,o — Lbulge relation 
is probably an overestimate, the total BH mass density 
from the M.,o — Lbulge and the galaxy LF may be actually 
not much different from that estimated from the M #i o — a 
relation and the galaxy velocity dispersion distribution 
function. 

3. THE QSO/AGN LF IN THE OPTICAL AND HARD X-RAY 

BANDS 

3.1. The optical QSO LF 

The optical QSO LF was first estimated by Schmidt 
(1968) and Schmidt & Green (1983), and it has been 
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Fig. 1. — The local BHMF obtained from the velocity-dispersion/luminosity distribution function of nearby galaxies and the MBH mass 
versus velocity-dispersion/luminosity relation (see details in § 2). Upper left panel: the solid lines represent MmfiTiM, [M m fi, to) obtained 
from the Af.o — cr relation given in Lauer et al. (2007a, here Lauer07), with assumed intrinsic scatters Aj^f. — a = 0.3, 0.2, and dex 
from top to bottom, respectively; while the dotted lines represent the M m ,onM m (M» o,to ) obtained from the M.,o — cr relation given in 
Bernardi et al. (2007, here Bernardi07) with the same assumed intrinsic scatters. The velocity-dispersion distribution is obtained from 
equation (8). Upper right panel: BHMF estimated from the galaxy luminosity function and the M.,o ~ ^bulgo relation with assumed 
intrinsic scatter Am. — z, bulgo = 0.4, 0.3 and dex from top to bottom, respectively. The solid lines represent the BHMF obtained from 
the M. o — ^bulgo relation given in Lauer07, while dotted lines for the relation given in Bernardi07. The bulge luminosity of a galaxy is 
used here by modifying the galaxy luminosity function to the bulge luminosity function as shown in § 2.2. Bottom left panel: comparison 
of the BHMF obtained from the M. o — cr relation and that obtained from the M. o — i-bulge relation for the relations estimated either 



in Lauer07a or in Bernardi07. The intrinsic scatter for the M.,q — a (or Af. o — ^bulge) relation is assumed to be Ajf ( . 



0.3 dex (or 



0.4 dex), which is taken as the most probable number in this paper. Bottom right panel: similar to the bottom left panel 



but with corrections for bias as suggested by Tundo et al. (2007, here Tundo07). 



investigated extensively since then. The shape and evo- 
lution of the QSO LF has been well, though not perfectly, 
constrained due to recent surveys with unprecedentedly 
large redshift and luminosity spans (e.g., Boyle et al. 
2000; Wolf et al. 2003; Croom et al. 2004; Richards et 
al. 2005, 2006a; Jiang et al. 2006; Fontanot et al. 2007; 
Siana et al. 2007). Using a sample of more than 15,000 
QSOs at redshift z < 2.5 from 2Qz and 6Qz, Croom et 
al. (2004) obtained the binned QSO LF ^/ M (M h]th Zj) 
over the range 0.4 < z < 2.1 and the magnitude range 
Mb, < —22.5, where Mbj,i is the ith bin of the abso- 
lute magnitude and Zj is the jth bin of the redshift. For 
some high-redshift bins, the binned QSO LF at low lu- 



minosity is not available because of the flux limit of the 
surveys. The time integral of the QSO LF can be es- 
timated through direct summation by multiplying the 
binned QSO LF by the cosmic time duration as 



T' 



E 

3 



*M(M bj ,i,^) A *(^)) 



(13) 



where At(zj) is the cosmic time interval corresponding 
to the redshift bin Zj, is assumed to be outside 
observational bins, and the prime ' indicates the value 
obtained by summation over bins — in contrast the vari- 
able T without prime (see eq. 15) represents the time 
integral of a continuous fit to the QSO LF. These sum- 
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mations only give lower limits to the time integral of 
the QSO LF because the binned QSO LF, especially in 
the low-luminosity bins, does not extend to high enough 
rcdshift to include all QSOs. Richards et al. (2006a) ob- 
tained the binned QSO LF over a larger redshift range 
(0.3 < z < 5) using a homogeneous statistical sample of 
15,343 QSOs drawn from SDSS Data Release 3. Unfor- 
tunately the SDSS survey is shallow so the binned QSO 
LF can only be determined at the bright end. As a com- 
plement to the above estimates, the QSO LF for faint 
QSOs over the range 1.2 < z < 4.8 was estimated by 
Wolf et al. (2003) using the COMBO-17 data; by Jiang 
et al. (2006) over the range 0.5 < z < 3.6 using a deep 
survey of faint QSOs in the SDSS; by Fontanot et al. 
(2007) in the redshift range 3.5 < z < 5.2 by combin- 
ing the data from the Great Observatories Origins Deep 
Survey (GOODS) and the SDSS; and by Siana et al. 
(2007) at the redshift range 2.83 < z < 3.44 using the 
data from the Spitzer Wide-area Infrared Extragalactic 
(SWIRE) Legacy Survey. In Figure 2, the direct sum- 
mations (eq. 13) are shown for the binned QSO LF from 
Croom et al. (2004, blue triangles), Wolf et al. (2003, ma- 
genta circles) and Richards et al. (2006a, green squares), 
respectively (the Mi magnitude in Richards et al. 2006a 
and M145 magnitude in Wolf et al. 2003 are all con- 
verted to M B magnitude by M B ~ M;(z = 2) + 0.80 and 
M B = M145 + 1.75, see Richards et al. 2006a and Wolf 
et al. 2003). At the high- luminosity end, the estimate 
from Croom et al. (2004) is substantially smaller than 
that from Richards et al. (2006a) which emphasizes the 
significance of the contribution from high-redshift QSOs. 
At the low-luminosity end, the estimates from Richards 
et al. (2006a) are smaller than those from others because 
the Richards et al. (2006a) sample is shallower and the 
majority of faint QSOs are not included. 

We combine these binned QSO LFs obtained by differ- 
ent surveys over different redshift and luminosity ranges 
(Croom et al. 2004; Wolf et al. 2003; Richards et al. 
2006a; Jiang et al. 2006; Fontanot et al. 2007; Siana et 
al. 2007), to cover luminosity and redshift ranges as large 
as possible. The basic rule is that the binned QSO LF 
from the largest sample are adopted at each redshift bin 
with data available and interpolations of the data points 
over magnitudes at a given redshift are used. The red 
points in Figure 3 are the estimated T^ b qsq with mean 
magnitude corresponding to that in Croom et al. (2004). 
At bright magnitudes, most of the points cover the range 
0.4 < z < 5 but the two points with faintest magnitudes 
only cover the range 0.4 < z < 1.0. In addition, the 
five green squares represent the brightest QSOs obtained 
from Richards et al. (2006a) only and are consistent with 
the trend of the red points. 

The optical QSO LF is frequently fitted with a double 
power law: 

$> M (M,z) = 100 4(A+1)[M _ M , (z)] - 10 o.4(/3 2 + l)[M-M*( z )] ' 

where VP m(M, z)dM is the comoving number density of 
QSOs with absolute magnitude in the range [M, M+dM] 
at redshift z. That is, the evolution of the QSO LF can 
be characterized by three functions of redshift: the slopes 
at both the high-luminosity (/3i) and the low-luminosity 
ends (f3 2 ) and the break luminosity (corresponding to 



M*). Boyle et al. (2000), Croom et al. (2004), Richards 
et al. (2005), and Jiang et al. (2006) all use this func- 
tional form to fit their data sets from 2dF and SDSS, 
except that Jiang et al. (2006) introduced additional den- 
sity evolution to the QSO LF at high redshift (z > 2.0). 
Adopting their best-fit models, the time integral of the 
QSO LF, 

/dt 
* M (M,z)—dz, (15) 

is obtained by integrating the QSO LFs over the range 
< z < 8. This function is shown in Figures 2 and 
3. There are some differences in the model parameters 
among the best-fit models for different samples. For ex- 
ample, Croom et al. (2004) obtained a slope of (3 2 ~ 1-09 
at the faint end (blue line), but Richards et al. (2005) 
obtained a steeper slope (f3 2 ~ 1.45; green line) us- 
ing a sample from the 2dF-SDSS LRG and QSO sur- 
vey (2SLAQ) with a flux limit of one magnitude fainter, 
which is roughly consistent with that obtained by Boyle 
et al. (2000) (/3 2 ~ 1.58; red line). Jiang et al. (2006) 
also obtained a shallower slope (fl 2 ~ 1-25; cyan line) 
with a deep survey in the SDSS, which is similar to 
that (/3 2 = 1.24) found by Hunt et al. (2004) at redshift 
z ~ 3. At high rcdshift (z > 3), the estimate of the faint- 
end slope by Fontanot et al. (2007) is consistent with 
/?2 = 1-45 but may have a high probability to be as steep 
as (3 2 = 1.71, and Siana et al. (2007) obtained (3 2 = 1.42, 
which is not inconsistent with values measured at lower 
redshift (e.g., Richards et al. 2005; Boyle et al. 2000). 
The differences in (3 2 are the primary reason for the dif- 
ferences in Tmb.qso at the faint end (see Figs. 2 and 3). 
(Below we choose (3 2 <~ 1.45 as the best estimate of the 
faint end of the QSO LF in § 5.) At the high luminosity 
end, the direct summations from the combination of the 
binned QSO LFs (according to eq. 13), which should be 
a lower limit to the time integrals, are quite consistent 
with the integration obtained from extrapolations of the 
best-fit analytic models, which may suggest that the es- 
timates of 7m b ,qso' a t l eas t at the high-luminosity end, 
are quite secure. 

3.2. X-ray AGN LF 

The advantage of counting QSOs/AGNs in X-rays is 
that relatively low-luminosity AGNs and obscured (type- 
2) AGNs, which may be missed in optical surveys, can 
be unambiguously detected in deep X-ray surveys even 
at large redshift. Although the number of QSOs/AGNs 
observed in X-rays (< 1000) is still substantially smaller 
than that observed in the optical band (> 10 4 ), the X- 
ray AGN (XAGN) LF can be estimated with consider- 
able accuracy (e.g., Ueda et al. 2003; La Franca et al. 
2005; Hasingcr et al. 2005; Barger et al. 2005; Silverman 
et al. 2008). Ueda et al. (2003) estimated the hard X-ray 
(2-10 kcV) LF (HXLF), which is assumed to repre- 
sent the total X-ray LF of unobscured plus Compton- 
thin AGNs, from a complete sample with <~ 257 sources 
observed by ASCA (but most of their sources have red- 
shift z < 3). La Franca et al. (2005) estimated the HXLF 
using a combined sample with 508 sources with redshift 
z < 2.5. With the data from Chandra deep surveys, 
Barger et al. (2005) extended the estimate of the HXLF 
(2-8 kcV) to higher redshift (3 < z < 5) but with 
large uncertainties at this redshift range. Combining 
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the published data from deep surveys by Chandra (i.e., 
CDF-North, CDF-South) and XMM-Newton (Lockman 
Hole) and rare luminous sources from the Chandra Mul- 
tiwavelength Project, Silverman et al. (2008) estimated 
the HXLF (2 - 8 kcV) at rcdshift 3 < z < 5 with much 
smaller uncertainties. The soft X-ray (0.5 — 2 keV) LF 
recently computed by Hasinger et al. (2005) is assumed 
to represent the unobscured type-1 AGNs. Gilli et al. 
(2007) demonstrated that the soft X-ray LF obtained 
by Hasinger et al. (2005) is actually consistent with the 
HXLF obtained by Ueda et al. (2003) and La Franca et 
al. (2005) by assuming a distribution of absorption col- 
umn densities. However, the bolometric correction (BC) 
for the soft X-ray band is much more uncertain than that 
in the hard X-ray band, so we shall not consider the soft 
X-ray LF further in this paper. The shape and evolution 
of the X-ray LF in both hard X-ray and soft X-ray bands 
can be described by the "luminosity-dependent density 
evolution" model (e.g., Ueda et al. 2003; La Franca et al. 
2005; Silverman et al. 2008; Miyaji et al. 2000; Hasinger 
et al. 2005): 

d^(< L x ,z) #(<L^z = 0) 
*i og L(logi x ,z) = JT7r7 e(z), 



d log L x 



d log L x 



(16) 

where \t , i og L(logLx, z)d\ogLx is the comoving number 
density of QSOs with logarithm X-ray luminosity in the 
range [logLx, logL x +d\ogL x ], 



diP(<L x ,z = 0) 
d log L x 



A 



Lx 



e(z) = 



_((1 + z)p\ 



e{z c )[(l + z)/{l + z c )Y», 



Lx 

z < z c 
z> z c . 



and 



3 (L 



x 



z* c (Lx/L a y 



Lx > L a , 
Lx < L a . 



(17) 
(18) 

(19) 



To estimate the time integrals of the HXLF, we will use 
the HXLF obtained by La Franca et al. (2005) as their 
AGN sample is larger than that in Ueda et al. (2003) 
and that obtained in Silverman et al. (2008) as their 
X-ray LF extends to redshift z <~ 5. In Figure 4, the 
direct summations obtained by multiplying the binned 
HXLF by the cosmic time duration in each luminosity 
bin with available data in each redshift bin are shown 
as green and red points for the HXLFs obtained in La 
Franca et al. (2005) and Silverman et al. (2008), respec- 
tively. The 2—8 kcV luminosity in Silverman et al. (2008) 
is converted to the 2 — 10 keV luminosity by assuming 
a photon index of 1.9. The time integral obtained by 
integrating the HXLF over redshift < z < 8 (with 
extrapolation of the HXLF to high redshifts and high lu- 
minosities) is shown in Figure 4 by adopting the best-fit 
"luminosity-dependent density evolution" model of the 
HXLF in Ueda et al. (2003) (blue line), La Franca et al. 
(2005, model 4 in table 2) (red line), and Silverman et al. 
(2008, model C in table 4; green line), respectively. In 
Figure 4, the direct summations obtained by multiplying 
the binned HXLF by the cosmic time duration, represent- 
ing the lower-limits to the time integrals of the HXLF, 
are quite consistent with the time integrals obtained by 
integrating the best-fit X-ray LF models, which might 
suggest that the majority of X-ray AGNs have been cov- 
ered by current observations although the HXLF from 



La Franca et al. (2005) does not cover redshift z > 2.5 
and the sample of Silverman et al. (2008) lacks high- 
luminosity AGNs. At the low- luminosity end, the time 
integrals obtained from the Silverman et al. (2008) HXLF 
is smaller than that from La Franca et al. (2005) by a 
factor of ~ 2, which may be due to the selection bias of 
the magnitude limits in the survey of the Silverman et 
al. (2008) sample. Hereafter we take the estimates ob- 
tained from La Franca et al. (2005) at the low-luminosity 
end (Lx < 10 43 5 erg s" 1 ) as the best estimates, while at 
middle and high luminosities both the estimates from La 
Franca et al. (2005) and Silverman et al. (2008) are taken 
into account. 

The X-ray cosmic background at a few to 100 keV is 
believed to be produced by the integrated emission from 
AGNs (e.g., Comastri et al. 1995). Using the synthesis 
model to reproduce the observed X-ray background, a 
population of Compton-thick AGNs is required to match 
the high energy (at <~ 30 — 40 keV) X-ray background 
spectrum as measured by HEAO-1 (e.g., Gilli et al. 2007). 
The number density of these Compton-thick AGNs is es- 
timated to be at most ~ 30% of the total population at 
Lx > 10 43 5 erg s -1 and not larger than 45% at lower lu- 
minosity (e.g., Gilli et al. 2007; Miiller & Hasinger 2007). 
A low-limit of the fraction of Compton-thick AGNs to the 
total population is probably ~ 10 — 15%, which is set by 
the current observations by INTEGRAL and Swift for 
bright AGNs (Markwardt et al. 2005; Beckmann et al. 
2006); and locally the fraction of Compton-thick AGN is 
found to be less than 20% by Sazonov et al. (2007). Cur- 
rent observations of the Compton-thick AGN population 
are insufficient to give its (luminosity) distribution func- 
tion. We will discuss the contribution of this population 
to the time-integral of AGN LF and its effect on model 
parameter, but do not go into details of the Compton- 
thick population in the models in § 4. 



3.3. The BC in the optical and hard X-ray bands 

The BC of a QSO is usually defined by C v = 
L^oxjiyLv), where vL v is the energy radiated at the 
central frequency v of a specific band. Based on obser- 
vations from optical to hard X-rays, Elvis et al. (1994) 
constructed the spectral energy distributions (SEDs) for 
several tens of QSOs and estimated the BC in the B 
band, which is about 11.8 ± 4.3. Considering that the 
infrared bump in the Elvis et al.'s SED templates was 
probably due to reprocessing of UV to X-ray photons by 
the dusty torus rather than the intrinsic emission from 
the central nuclei, Marconi et al. (2004) obtained that 
the BC at the B band is 7.9 ± 2.9. Based mainly on an 
anti-correlation between the optical-to-X-ray spectral in- 
dex (a ox ) and the 2500A luminosity (e.g., Vignali et al. 
2003; Strateva et al. 2005; Steffen et al. 2006), Marconi 
et al. (2004) and Hopkins et al. (2007a) re-calibrated the 
SED and argued that the BC is luminosity-dependent. 
The BCs were derived by Marconi et al. (2004) as 

log[L bo i/L 2 _ 10 kcV ] = 1.54 + 0.24L + 0.012L 2 -0.0015L 3 , 

(20) 

\og[L hoX /v B L^] = 0.80 - 0.067L + 0.017L 2 - 0.0023L 3 , 

(21) 
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Fig. 2. — The time integral of the QSO luminosity function in 
the optical band (B). The solid triangles (blue), solid circles (ma- 
genta) and solid squares (green) represent the direct summation of 
the binned QSO LF (eq. 13) obtained from Croom et al. (2004) 
over the redshift range 0.35 < z < 2.1, Wolf et al. (2003) over the 
redshift range 1.2 < z < 4.8 and Richards et al. (2006a) over the 
redshift range 0.3 < z < 5.0, respectively. High-luminosity QSOs 
are under-represented in Croom et al. (2004) because of the redshift 
limit (z < 2.1) and low-luminosity QSOs are under-represented in 
the sample in Richards et al. (2006a) because of incompleteness at 
the faint end (especially at high redshift). The blue, cyan, red, and 
green lines represent the time integrals obtained from the fitting 
formulae in Croom et al. (2004), Jiang et al. (2006), Boyle et al. 
(2000), and Richards et al. (2005), respectively. The time inte- 
grals obtained from different fitting formulae are consistent at the 
high-luminosity end but show substantial discrepancies at the low- 
luminosity end (Mb > —23) mainly because of the uncertainties 
in the faint-end slope of the QSO LF. For example, the faint-end 
slope is estimated to be ~ 1.09 by Croom et al. (2004), ~ 1.25 in 
Jiang et al. (2006), but ~ 1.58 by Boyle et al. (2000) and ~ 1.45 
by Richards et al. (2005). See also in Fig. 3. 



where L = log -Lboi — 12 and Lboi is the bolometric lumi- 
nosity in units of L Q . Hopkins et al. (2007a) found 

fei /j \ fea 

n|— I Wt^M , (22) 



^bol 



-'band 



10 10 L 



10i°L Q 



with (ci, kx,c 2 ,k 2 ) given by (6.25, -0.37, 9.00, -0.012) for 
iband = Lb and (10.83, 0.28, 6.08, -0.020) for L band = 
L2—10 koV- The scatter in BCs given by equation (22) is 



^(ibd/i^d) 



= ai(i bol /10 9 L c . 



02, 



(23) 



where (cti, (3, cr 2 )=(0.08, -0.25, 0.06) in the B band and 
(0.06, 0.10, 0.08) in the hard X-ray. The BC in hard 
X-ray given by Hopkins et al. (2007a) is 30% larger than 
that given by Marconi et al. (2004), and the BC in the 
B band given by Marconi et al. (2004) is smaller than 
that given by Hopkins et al. (2007a) by a factor of 1.5 
(or 1.8) at Lboi = 10 10 L© (or L bo i = 10 14 i Q ). In this 
paper, we adopt the BCs for the X-ray and B bands and 
associated scatters obtained by Hopkins et al. (2007a). If 
the BCs given by Marconi et al. (2004) were adopted, the 
efficiency e should be systematically smaller than that 
obtained below in § 4 by a factor of ~ 1.3 in order to 
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Fig. 3. — The time integral of the QSO LF in the optical band 
(-B). Similar to Fig. 2, but the points are obtained by combin- 
ing the binned optical LFs given by different surveys over different 
redshift and luminosity ranges. At each redshift bin with data 
available, the binned optical LFs obtained from the largest sample 
are adopted, and interpolations of the data points over magni- 
tudes at a given redshift are used. The points, which should be 
lower limits to the time integrals of the QSO LF, are quite con- 
sistent with that obtained from the fitting formulae of the QSO 
LF at Mb < —24 (solid lines). At the faint end, the direct sum- 
mations are substantially smaller than those estimated from the 
continuous fitting formulae with extrapolations to higher redshift 
and lower luminosities (which may be due to the incompleteness of 
the samples). The five green squares represent the brightest QSOs 
obtained from Richards et al. (2006a) only and they are consistent 
with the trend of the red points. The blue, cyan, red and green 
lines represent the time-integrals obtained from the fitting formu- 
lae in Croom et al. (2004), Jiang et al. (2006), Boyle et al. (2000) 
and Richards et al. (2005), respectively. 



match the time- integral of QSO/AGN LF obtained from 
observations with that inferred from the local BHMF. 

We note that Vasudevan & Fabian (2007) recently in- 
vestigated the SEDs of 54 AGN and found significant 
spreads in the BCs. Their results suggest a relationship 
between BCs in the X-ray band and Eddington ratios 
(see definition in § 4) in AGNs, with a transition at an 
Eddington ratio of ~ 0.1, below which the BC is typi- 
cally 15 — 25 for the 2 — 10 keV luminosities and above 
which the BC is typically 40 — 70. Their estimates of 
the BC for the optical band is approximately indepen- 
dent of Eddington ratio and roughly consistent with that 
obtained by Hopkins et al. (2007a). We also note that 
simple theoretical expectations of the BCs would be that 
it is not only the functions of Eddington ratios but also 
the functions of MBH masses because the SED of the 
disk emission depends on the MBH mass and Edding- 
ton ratio. In addition, the QSO/AGN variability in the 
hard X-ray is substantial while it is not significant in the 
optical band. The X-ray variability, typically a factor 
of ~ 1.5, introduces an additional scatter of ~ 0.13 dex 
to the BC for the hard X-ray band (see Tab. 2 in Va- 
sudevan & Fabian 2007). Since a quantitative relation 
between the BCs and the Eddington ratio is still prema- 
ture, we shall not consider the BCs as functions of the 
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Fig. 4.— The time integral of the X-ray AGN LF. The red solid 
circles represent the direct summations of the binned X-ray LF 
for the AGN sample over redshift range < z < 2.5 collected in 
La Franca et al. (2005), while the green solid circles represent the 
direct summations of the binned X-ray LF obtained from the AGN 
sample over redshift range 0.2 < z < 5.5 by surveys such as the 
Chandra Deep Field (CDF) described in Silverman et al. (2008, 
see details in § 3). For the low-luminosity points, the sample may 
be incomplete at high redshifts and thus those points may be only 
lower limits, especially for the sample in Silverman et al. (2008) 
which may suffer from a bias due to the optical magnitude limits 
in the survey. The red, blue, and green lines represent the time 
integrals of the X-ray LFs in the redshift range < z < 8 according 
to the fitting formulae obtained by La Franca et al. (2005), Ueda 
et al. (2003), and Silverman et al. (2008), respectively. At the 
high-luminosity end (Lx > 10 44 erg s — 1 ), the time integrals are 
quite consistent with the binned data and with each other, but at 
the low-luminosity end (Lx < 10 4s erg s — 1 ), the one obtained from 
Silverman et al. (2008) is smaller than that from La Franca et al. 
(2005) by a factor of > 2. The estimates obtained from La Franca 
et al. (2005) at the low-luminosity end (Lx ^ 10 43 5 erg s _1 ) is 
taken as the best one, while at middle and high luminosities both 
the estimates from La Franca et al. (2005) and Silverman et al. 
(2008) are taken into account. Compton-thick sources, which are 
hard to be observed even in the X-ray, are not included (for the 
contribution of the Compton-thick AGNs to the time-integral of 
X-ray luminosity function see discussions in § 3.2). 

Eddington ratio in this work but simply adopt equations 
(22) and (23) and include an additional scatter due to 
the X-ray variability. 

The time integral of the QSO LF at any given wave- 
band Y can also be inferred from the local BHMF as 
follows, provided that the BC at this band is known, 

/•oo /*oo 

T™x d AGN= / dL ho} / n M .(M.,o,to)rit(M.,o) 
Jo Jo 

P(L hol \M., )P(L Y \L hol )dM., , (24) 

where Ly is the luminosity at the Y-band, and 
P(Ly \Lho\) is the probability distribution of Y-band lu- 
minosity for QSOs/AGNs with bolometric luminosity 
Lboi and is determined by the BCs and their scatters. 

4. SIMPLE MODELS FOR THE LUMINOSITY EVOLUTION 
OF INDIVIDUAL QSOS 

In this section, we introduce three simple models 
for the luminosity/accretion rate evolution of individual 



QSOs. These models are assumed to represent the lumi- 
nosity/accretion rate evolution averaged over an inter- 
mediate timescale substantially smaller than the lifetime 
of individual QSOs, but much longer than certain de- 
tails of the evolution such as the short time variation, 
etc. The parameters involved in these models will then 
be constrained by observations of the local BHMF and 
the QSO/ AGN LF through the extended Soltan argu- 
ment (eq. 1). Because X-ray surveys are more complete 
than optical surveys in the sense that obscured AGN can 
be detected in X-ray surveys, we will compare the time 
integrals obtained from the X-ray LF with that inferred 
from the local BHMF in this section, and then use the 
time integral of the optical QSO LF to give constraints 
on obscured AGN fraction in the optical band in § 5. 

4.1. Several fiducial parameters 

We first summarize several fiducial parameters in- 
volved in the models below. 

• The "Eddington luminosity" is a characteristic lu- 
minosity at which radiation pressure on free elec- 
trons balances gravity: 



^Edd 



(M.) = 



4nGM.m p c 



or 



1.26 x 10 



40 



/ M. 



Vio 8 m 



erg s 



(25) 



where G is the gravitational constant, m p is the 
mass of a proton, and ctt is the cross-section of 
Thompson scattering. The Eddington luminosity is 
frequently assumed to be the maximum luminosity 
of any object of mass M m . 

• Corresponding to the Eddington luminosity, the 
"Eddington accretion rate" is defined by: 



M°° 

iw acc,Edd- 



L 



Edd 



M Q yr" 



(26) 



where e is the mass-to-energy conversion efficiency; 
and the Eddington growth rate of a MBH is 



M. 



Edd 



(1 - e)M ( 



acc.Edd" 



(27) 



The efficiency e is predicted to be in the range 
~ 0.04 — 0.31 in the thin disk accretion models, 
depending on the spin of the MBH [e = 0.057 
for a Schwarzschild BH, and 0.31 (0.04) for a pro- 
grade (retrograde) rotating accretion disk around 
a Kerr BH with the dimensionless spin parameter 
a ~ 0.998, the upper limit of BH spin if the BH is 
spun up by accretion; Thorne 1974]. Currently, the 
spin of MBHs is difficult to measure directly. The- 
oretical studies of the spin evolution of MBHs show 
that MBH spin may reach an equilibrium point for 
most of its lifetime considering both accretion and 
merger processes (e.g., Lu et al. 1996; Gammie et 
al. 2004; Shapiro 2005; Volonteri et al. 2005; Haw- 
ley et al. 2007; Noble et al. 2008; Hughes & Bland- 
ford 2003). This equilibrium value is - 0.7 - 0.9 
and corresponds to an efficiency e ~ 0.10 — 0.20 
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(e.g., Gammie et al. 2004; Shapiro 2005; Hawley et 
al. 2007). If the accretion rate of a MBH is less 
than the Eddington rate by a factor much larger 
than 100 (e.g., to = M./M. Edd < 10~ 3 ), the MBH 
may accrete material via the Advection Dominated 
Accretion Flow (ADAF) with very low efficiency, 
e « 0.1 (e.g., Narayan & Yi 1994), or via a mode 
described by the Advection Dominated Inflow and 
Outflow scenario (ADIOS, Blandford & Begelman 
1999) with most of the accretion material blown 
away. The contribution from these very low ef- 
ficiency modes to the observational range of the 
time integral of the QSO / AGN LF is negligible and 
MBH growth may also be very inefficient in this 
low-accretion rate mode. In this paper, we will not 
consider this complication but assume that e is a 
constant that is neither directly nor indirectly re- 
lated to the BH mass M. and the accretion rate, 
as e is probably mainly determined by the spin of 
the central BH in the thin-disk accretion mode. (A 
more detailed study of the growth of MBHs should 
simultaneously consider the spin and mass evolu- 
tion of MBHs.) 

• If a MBH-disk accretion system accretes material 
via the Eddington accretion rate and radiates with 
luminosity Lboi, the mass of the MBH is 



^•,Edd(£bol) = 



1.26 x 10 38 erg s" 



Mr. 



(28) 



• The Salpeter timescale is defined as the time for 
a MBH radiating at the Eddington luminosity to 
c-fold in its mass: 



M. 



M. 



= 4.5 x 10 7 



,Edd 



0.1(1 -e) 



yr. 



(29) 



If the accretion rate is only a fraction A of the 
Eddington accretion rate, then the timescale for 



a MBH to e-fold its mass is Tg 



4.2. Model (a) 

The mass of MBHs in QSOs may be estimated by using 
the virial mass estimator(s), i.e., using the width of broad 
emission lines and the empirical relation between the op- 
tical luminosities and the sizes of broad line regions esti- 
mated from reverberation mapping studies (e.g., Wandel, 
Peterson & Malkan 1999; Kaspi et al. 2000; Vestergaard 
2002; Kaspi et al. 2005; see also discussion of uncertain- 
ties, e.g., in Krolik 2001), and hence the Eddington ra- 
tio may be estimated (e.g., Woo & Urry 2002). Recent 
studies by Kollmeier et al. (2006) on a sample of QSOs 
using the virial mass estimator(s) have suggested that 
the Eddington ratios (to = Lboi/^Edd) in QSOs, may be 
consistent with a single value, and the best estimates of 
the mean value of to is around 10~ 06 for all redshifts 
and luminosities. Using a large sample of QSOs from 
SDSS, Shen et al. (2007) investigate the Eddington ra- 
tio distribution in QSOs over a range of redshifts and 
luminosities, however, their results show that the mean 
value of the Eddington ratio is a function of redshift and 
luminosity and it ranges from 10 -1 - 1 and lO" 06 . Netzcr 
et al. (2007) also argue that the to distribution is not 
consistent with a single value and the conclusion that 



a single to applies to all QSOs/AGNs might be due to 
some unknown selection effects. Ignoring this concern, 
for the moment, we assume that all MBHs in QSOs ac- 
crete material at a constant normalized rate to = A, i.e., 
M. = AM. Edd while the QSO is "on". The luminosity 
evolution is 



L ho i(M. fi , t) = XL Edd (M, i0 ) exp 



Tlt(M. ;0 ) 



for < t < r lt (M. ;0 ). 



(30) 



This model involves three parameters (e, A, £), where 
£ = Tit(M., )/Tg = ATit(M. !0 )/Ts; and these three pa- 
rameters solely determine the growth history of individ- 
ual MBHs. For MBHs with present-day mass M,.o, the 
probability distribution of the nuclear luminosity in their 
evolutionary history (eq. 3) is 



P(L bol |M. j0 ) 



fgj_ 



where 



f Q 



if M -^ (Lb ° l} < M., < 
otherwise, 



bol 



M. 



(31) 



l(ibol) 



exp(0, 



and the present-day mass of a MBH is related to its initial 
mass M, t i at the time of nuclear activity being triggered 
by M.,o = M #)i exp(£). 

For a given set of parameters (e, A, £), we calculate 
the time integrals of the XAGN LF, T L ™ ^ AGN (or 
^Mb"xagn ) > using equation (24). To do this, the local 
BHMF is chosen to be the one estimated by using the 
M #i o — tr relation by Lauer et al. (2007a) as the refer- 
ence BHMF in this paper (see the solid blue line in the 
right bottom panel of Fig 1). For given BHMF, BCs, 
and A, the normalization of the inferred time integrals of 
XAGN LF, i.e., ^"xagn > is proportional to e/(l - e) 
through t' s , which can vary by a factor of 10 for the typ- 
ical range of e, 0.04-0.31. If £ is substantially smaller 
than 1, T^xXAGN is also proportional to £ because the 
range of the integration limits over the BH mass is quite 
small and thus approximately proportional to £, and the 
shape of the inferred time integrals of the QSO LF is de- 
termined by the shape of the local BHMF. In this case, 
there is some degeneracy between the parameters £ and 
e if £ < 1. However, this degeneracy does not exist if £ is 
substantially larger than 1 (i.e., if the growth of MBHs is 
dominated by accretion processes) [which is also true for 
models (b) and (c) below] , as xagn ^ s insensitive to £ 
at the high- luminosity end and increases only slowly with 
increasing £ at the low-luminosity end. For example, the 
predicted T^xagn f° r tnc case of £ = 2 (but fixed e 
and A) at the low- luminosity end (Lx < 10 43 crg s _1 ) is 
smaller than that for the case of £ = 10 (with the same e 
and A) by a factor of ~ 1.2 — 1.3, and T^xagn for these 
two cases are almost the same at the high-luminosity end 
(L x > lO^ergs" 1 ). 

As shown in Figure 5, A should be in the range from 
0.5 to 1 in order to match the observations at high lu- 
minosity (L x > 10 44 erg s" 1 ) with T£™^ AGN , while it 
should be close to 0.1 in order to match the observations 
at lower luminosities (Lx < 10 43,5 erg s _1 ). According 
to Figure 5, we conclude that the inferred time integrals 
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of the XAGN LF cannot match the observations at both 
the low-luminosity and the high-luminosity ends simul- 
taneously if all MBHs accrete material at a single rh = A. 

We should note here that T^XAGN may well fit the 
observations if A is arbitrarily assumed to be an increas- 
ing function of M #) o (cf., the Eddington ratio may be 
redshift-dependent and thus mass-dependent since sta- 
tistically MBHs with larger M. : o formed earlier, see 
Shankar et al. 2004, 2007). However, the assumption that 
all low-mass MBHs need to accrete material via lower 
Eddington ratios may be not realistic/physical because 
(1) some low-mass MBHs, such as the one in NGC 3079 
(with MBH mass - 2x 1O 6 M ) or NGC 1068 (with MBH 
mass ~ 8 x 1O 6 M ), do accrete material with a rate close 
to the Eddington limit and have massive accretion disks 
with mass comparable to the MBH mass (Kondratko et 
al. 2005; Lodato & Bcrtin 2003); and (2) there is no clear 
physical reason for the low-mass MBHs to accrete mate- 
rial via smaller Eddington ratios compared to high-mass 
MBHs if they also obtained their mass mainly from ac- 
cretion. Therefore, we do not pursue the possibility that 
the Eddington ratio is constant for each AGN with the 
same MBH mass but an increasing function of the MBH 
mass. 

4.3. Model (b) 

A more realistic model would be that the growth of 
MBHs involves two phases after the nuclear activity is 
triggered (see the discussions in Small & Blandford 1992, 
Blandford 2003, and Yu & Lu 2004a). In the first phase, 
there is plenty of material to supply the MBH growth; 
however, MBHs may not be able to accrete as fast as ma- 
terial fueling allows because the accretion process may 
be self-regulated by the Eddington limit. With the de- 
cline of the material supply, the MBH growth enters the 
second phase and the nuclear luminosity in which the 
limiting factor is the fuel supply and accretion rate are 
expected to decline to below the Eddington limit. 

After the nuclear activity of a MBH is triggered at 
cosmic time we assume that the MBH accretes mate- 
rial via the Eddington accretion rate for a time-period of 
Tp = £ts, hence its mass increases to M #i p and its lumi- 
nosity approaches a peak of Lp(M #j0 ) = L Edd (M. tP ) at 
time tp = ti + Tp. The nuclear luminosity in this phase 
increases with time as 

Lboi(r) = £ Bd d(Af.,p) exp \ J T ^ F J > < t < t p , 

(33) 

where r = t — U is the age of the QSO since the nuclear 
activity was triggered. 

In the second phase, we assume that the evolution of 
the nuclear luminosity (or accretion rate) declines expo- 
nentially as (e.g., Haehnelt et al. 1998; Haiman & Loeb 
1998): 

£boi(M.,o,r) = 

I L Bdd (M. tP ) exp , for r P < r < t p + r}T D , 

{ 0, for t > Tp + t)Td, 

(34) 

where to = (ts is the characteristic decay timescale 
of the nuclear luminosity. We assume that QSOs be- 
come quiescent when the nuclear luminosity declines by 
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Fig. 5. — Comparison of the time integral of the X-ray AGN 
luminosity function (XAGN LF) and that inferred from the local 
BHMF by adopting model (a) in § 4, i.e., assuming that all MBHs 
accrete material via a fixed Eddington ratio (A). The symbols and 
the red line are the same as in Fig. 4. The black lines represent 
the inferred time integrals of the XAGN LF from the local BHMF 
with e = 0.14, £ = 10, and A = 1 (solid line), 0.5 (dotted line), 0.25 
(short-dashed line), and 0.08 (long-dashed line), respectively. The 
bottom panel shows the inferred time integrals of the XAGN LF 
from the BHMF compared to the prediction of the fitting-formula 
of XAGN LF obtained by La Franca et al. (2005). Note that the 
estimate of the highest-luminosity point at Lx = 10 46,5 erg s — 1 
cannot fit into any model, which might be partly because this point 
is estimated in La Franca et al. (2005) from only two AGN with 
luminosity ~ 10 461 erg s^ 1 in the bin 10 46 — 10 47 erg s — 1 . We use 
these two AGNs to give an estimate on the space density of AGN 
at ~ 10 46 1 erg s" 1 in a bin 10 46 — 10 46,2 erg s~ 1 and show it as the 
open circle in this figure. The estimates obtained from La Franca et 
al. (2005) at the low-luminosity end (Lx < 10 43 ' 5 erg s — 1 , the red 
points) are adopted as the best one since the low-luminosity data 
from Silverman et al. (2008) may suffer from selection bias, while at 
middle- and high-luminosities (10 43 5 erg s _1 < Lx < 10 46 erg s — 1 , 
the green points) both the estimates from La Franca et al. (2005) 
and Silverman et al. (2008) are taken into account. As shown in 
this Figure, the inferred time-integrals of XAGN LF cannot match 
the observations simultaneously at both the low-luminosity and 
high-luminosity ends. 

a factor of exp(— rf) compared to the peak luminosity 
LEdd(Mm.p), so there is a cutoff of the nuclear luminos- 
ity at Tn = Tp + r\Tu = (£ + ?;C) T s hr equation (34). 
The factor rj is set to — ln(10~ 3 ) = 6.9 here, since af- 
ter decreasing by a factor of 10~ 3 in accretion rate, the 
accretion mode may change from the efficient thin-disk 
accretion to the inefficient advection dominated accretion 
modes and the nuclear luminosity of MBHs even with a 
high mass ~ 10 9 M Q will become fainter than the lumi- 
nosity range (Mb < —20 or Lx < 10 42 erg s _1 ) of interest 
in this paper. With the assumption that all QSOs are 
quenched at present (i.e., to — ti — Tp ^> td), the MBH 
mass at the present day is 

M.,o =s (1 + C) M.,p = (1 + C) eM0 M ;i- (35) 



Precise Constraints on MBH Growth 



13 



This model reduces to model (a) with A = 1 if the second 
phase is not significant. 

For MBHs with present-day mass M, : q, the probability 
distribution of the nuclear bolometric luminosity in their 
evolutionary history (eq. 3) is 

P(L bol \M., )= fp l CfD -^-, (36) 

where 

fl if (l + C)M., Edd (L b oi) < M.,o < 
fp=\ (l + C)exp(OM. !Edd (L bol ), (37) 

t otherwise, 

and 

(I if (1 + C)M. Edd(iboi) < M.,o < 

f*>=\ 10 3 (l + C)A/ vEdd (L b oi), (38) 

1 otherwise. 

In model (b), we also have three parameters (e,£, C) to 
be constrained below. 

For any given set of parameters (e, £, £), we calcu- 
late the time integrals of the XAGN LF, 7^° X agn- Our 
calculations show that the dependence of T^xagn on 
parameters e or £ for a given £ is similar to that in 
model (a). For given e and £, a larger £ corresponds 
to smaller T^xagn at higher luminosities but larger 

T^xagn a ^ l° wer luminosities. As shown in Figure 6, 
£ must be around or smaller than 0.1 — 0.3 to match 
^Lx°xagn with observations at the high-luminosity end 
(Lx > 10 44 5 erg s" 1 ), but ( should be larger than 1 to 
match T^xagn with observations at the low-luminosity 
end (L x < 10 43 erg s" 1 ). It is unlikely that T^xagn 
inferred from any single £ with any fixed (e, £) can match 
observations simultaneously at both the high- and low- 
luminosity end. 

As discussed in model (a), £ and £ are not necessarily 
constants in model (b) but may be functions of M..o; or 
alternatively the ratio of the MBH mass M, t p at the peak 
luminosity to the final MBH mass M ti o may be a slowly 
increasing function of A/..o as proposed by Hopkins et 
al. (2006). The dependence of £ and £ on M, t o would be 
related to the assembly history of each MBH and the dis- 
tribution of seed BHs, which are poorly known. In model 
(b) , it is possible that 7/™° xagn can ma t cn observations 
at both the high- and low-luminosity ends if £ decreases 
with increasing M^n. But we will not go further to make 
this fit, for the same reasons given at the end of § 4.2. 

4.4. Model (c) 

The accretion rates in the second phase of QSOs, in 
which the luminosity decays, may be ultimately deter- 
mined by the evolution of the viscous accretion disk itself 
rather than galactic-scale dynamical disturbances. The 
disk accretion evolution may follow a self-similar solution 
(e.g., Pringle 1981; Lin & Pringle 1987; Cannizzo et al. 
1990; Pringle 1991), i.e., the accretion rate declines as a 
power-law of the QSO age (M, cx r -7 ), where the slope 
7 may be determined by the opacity law. The value of 7 
also depends on the binarity of the MBH surrounded by 
the accretion disk, for instance, 7 ~ 1.2 — 1.3 for the evo- 
lution of a disk around a single MBH (e.g., Cannizzo et al. 
1990), while 7 ~ 2.5 — 3.3 for a disk truncated by an outer 
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Fig. 6. — Similar to Fig. 5, but adopting model (b) in § 4, i.e., 
an initial accretion phase with rate set by the Eddington limit, 
followed by a phase with an exponentially declining accretion rate. 
The black lines represent the inferred time integrals of XAGN LF 
from the local BHMF with e = 0.16, £ = 10, and f = 0.1 (long- 
dashed line), 0.3 (short-dashed line), 1 (dotted line) and 3 (solid 
line), respectively. As shown in this Figure, the inferred time- 
integrals of XAGN LF cannot match the observations simultane- 
ously at both the low-luminosity and high-luminosity ends. 

secondary MBH (e.g., Lipunova & Shakura 2000). For 
the binary MBH system, however, the secondary MBH 
embedded in a disk surrounding the primary MBH may 
migrate inward and may merge with the primary MBH 
on a time-scale of 10 7 yr (e.g., Armitage & Natarajan 
2002; Escala et al. 2004, 2005), and thus the evolution 
of disk accretion associated with the binary MBH sys- 
tem (7 ~ 2.5 — 3.3) may not be sustained for a period 
substantially longer than 10 7 yr. The observed accretion 
rate distribution in local AGNs is found to be consistent 
with the self-similar evolution around a single MBH and 
7 ~ 1.26 ± 0.1 (Yu et al. 2005, see also King & Pringle 
2007b). In this paper we neglect the complications in the 
evolution of disk accretion due to possible binary MBHs, 
and assume 7 ~ 1.2 — 1.3. Below we introduce model (c) 
which is similar to model (b) but the nuclear luminosity 
in the second phase declines with time as a power law: 

iboi(M., ,r) = 

{-^Edd(-^»,P 
0, 



T + Tp-Tp 
TO 



, for Tp < T < Tp + T]T D , 
for T > Tp + TjTpt , 

(39) 

where td = C r s is the transition timescale from the first 
to the second phase. As in model (b), we assume that 
QSOs become quiescent when the nuclear luminosity de- 
clines by a factor of 10 3 compared to the peak luminosity 
LEdd(M, p) and afterwards the growth of MBHs is not 
significant. Thus r\ = 10 3/7 - 1 in this model. The MBH 



14 



Yu & Lu 



mass at a time r after the nuclear activity was triggered 
is 

' t — Tp s 



Ml = M #jP exp 



7"S 



(40) 



in the first phase, and is 



Ml = M. tP 



1 + 



c 



7-1 



td 



(41) 

in the declining phase. The present-day mass of a MBH 
is 

M. j0 ^X^.,P = Xexp(C)M., l , (42) 



1 + i-io^ 2^2 ^ T he slope can be 7 



(1 + 3.41C)M. iP , and (1 



where \ - , ^j- 

1.2, 1.3 and thus M #j0 ■ 
2.66C)M. ! p, respectively. 

In model (c), for MBHs with present-day mass M #j o 
the probability distribution of the nuclear bolometric lu- 
minosity in their evolutionary history (eq. 3) is 



P(Lboi|M., ) = 



JP + JD- { J 



1/7 



ibol 



where 



r 1 if xM vEdd ( J L bo i) < M.,o 

/p = S < xexp(£)M. iEdd (£ b oi) 

L otherwise, 



(43) 



(44) 



and 

/ D = / 1 if X^.,Edd(iboi) < M.,o < 10 3 xM vEdd (i bo i), 
1 otherwise. 

Besides the three parameters (e, £, £) involved in model 
(c), an additional parameter 7 is also involved, but 7 is 
fixed by assumption to be 1.2-1.3 here, if not otherwise 
specified, according to theoretical models on the long- 
term evolution of viscous disk (e.g., Pringle 1981; Lin 

6 Pringle 1987; Cannizzo et al. 1990; Pringle 1991) and 
recent observational constraints (e.g., Yu et al. 2005, see 
also King & Pringle 2007b). 

For any given £ in this model, the dependence of 
^LxXAGN 011 t ne parameters e and £ is similar to that 
in models (a) and (b). For given e and £, a larger £ is 
responsible for a smaller T^xagn at higher luminosi- 
ties but a larger T^xagn a * l° wer luminosities because 
a larger £ means that a larger fraction of the mass of 
MBHs is accreted via the second phase with Edding- 
ton ratios (substantially) smaller than 1. Although the 
growth history of both low-mass MBHs and high-mass 
MBHs is assumed the same in this model for fixed pa- 
rameters (r, £. (), apparently there are more objects with 
low Eddington ratios at the low-luminosity end but few 
objects with low Eddington ratios at the high-luminosity 
end. Detailed investigation of the Eddington ratio dis- 
tribution inferred from this model is discussed in § 6. 
As shown in Figure 7, T^xagn can match observations 
very well if C ~ 0.15 - 0.3 provided that e = 0.16, £ > 2 
and 7 ~ 1.2 — 1.3. With the parameter £ ~ 0.15 — 0.3 and 

7 ~ 1.2 — 1.3, the mass growth of MBHs at the accre- 
tion stage with Eddington ratio m < 1 ( or m < 0.1) 
is roughly a fraction ~ 0.2 — 0.5 (or <~ 0.1 — 0.3) of 



its final mass M. ; o, and this is compatible with the as- 
sumption that the disk mass is substantially less than 
the central MBH in the long-term evolution of disk ac- 
cretion model (e.g., Pringle 1981; Lin & Pringle 1987; 
Cannizzo et al. 1990; Pringle 1991); and MBHs obtained 
majority of their mass (> 80%) via a rate close to the 
Eddington limit (m > 0.1). With these parameters, we 

have 7i t (M.,o) ~ r s [£ + (10 3 ^ - 1)Q ~ (3 - 6) x 10 9 yr, 
and the period for MBH-accretion disk systems radiat- 
ing at luminosities larger than 10% of its peak luminos- 
ity, thus roughly rh > 0.1 (or m > 0.01), is only about 
(3 - 4)r s ~ (2 - 3) x 10 8 yr (or - 10 9 yr). Model (c) 
is based on detailed considerations of the evolution of 
disk accretion and appears to fit observations much bet- 
ter than models (a) and (b). Therefore, this model with 
three parameters (e, £, £) = (0.16, 10, 0.20) is set as the 
reference model in this paper. 

Considering of the uncertainty in the M #i0 — cr relation, 
the velocity dispersion distribution function and the time 
integral of XAGN LF, the error in the best-matched pa- 
rameter e is Se ~ 0.04. Note also that Compton-thick 
objects may be still missed in the hard X-ray surveys 
by La Franca et al. (2005). The fraction of Compton- 
thick objects should not be larger than 30% according 
to the X-ray background synthesis model (e.g., Miiller 
& Hasingcr 2007), and this would add additional uncer- 
tainty at most I 9 05 to e. To match the time integral of 
XAGN LF with the local BHMF, the efficiency e is re- 
quired to ~ 0.16 ± 0.04±°, 05 , and this range of e is fully 
consistent with theoretical expectations e ~ 0.10 — 0.20 
(e.g., Gammie et al. 2004; Shapiro 2005; Hawley et al. 
2007). The range of e (~ 0.12 - 0.25) constrained above 
corresponds to the spin parameter a in the range from 
~ 0.8 to 0.99 as the value of e is mainly determined by a 
with only an order of 10 — 20% or less uncertainty (e.g., 
Noble et al. 2008), which suggests that most MBHs in 
QSOs are indeed rapidly rotating Kerr BHs. It is worth 
to note that if we choose e = 0.21, the time integral 
of XAGN LF inferred from the local BHMF with pa- 
rameters (£ <~ 1, ( <~ 0.2) can still match the observa- 
tions well, but the time integral of XAGN LF at lumi- 
nosities Lx ^5 10 45 erg s _1 is overpredicted by ~ 40% 
if (£ <~ 10, C ~ 0.2) and the overpredicted part can 
be accounted for by the additional contribution from 
Compton-thick AGNs. Previous estimates of the effi- 
ciency include e > 0.15 (Elvis et al. 2002), e > 0.1 (Yu & 
Tremaine 2002; Yu & Lu 2004a), e - 0.04 - 0.16 (Mar- 
coni et al. 2004), e - 0.30 - 0.35 (Wang et al. 2006), and 
e - 0.06 - 0.11 (Shankar et al. 2004, 2007). 

The relatively high efficiency constrained above sug- 
gests that the majority of QSOs should not accrete ma- 
terial via the chaotic accretion scenario proposed by 
King & Pringle (2007a) to explain the rapid growth of 
MBHs in those QSOs at z > 6, in which MBHs spin 
down because of counter-alignments of their spin axes 
with accretion disk angular momenta and thus the effi- 
ciency reduces to a low value, close to the efficiency for 
Schwarzschild BHs. 

Note that there are some uncertainties in the intrinsic 
scatter in the M., — <J relation, which may mainly in- 
troduce some uncertainties to the parameter Q. A larger 
intrinsic scatter corresponds to more MBHs at the high- 
mass end and thus allows a larger £. But the uncertain- 
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Fig. 7. — Similar to Fig. 5, but adopting model (c) in § 4, i.e., 
an initial accretion phase with rates set by the Eddington limit, 
followed by a phase with power-law declining accretion rates, as 
might be set by self-similar long-term evolution of disk accretion. 
The black lines represent the inferred time integrals of XAGN LF 
from the local BHMF with e = 0.16, 5 = 10, and C = 0.15 (dotted 
line), 0.20 (solid line) and 0.30 (dashed line), respectively. This 
Figure shows that the time-integral of XAGN LF inferred from 
the local BHMF can well match observations simultaneously at 
both the low-luminosity and high-luminosity ends with suitable 
parameters. We take the model (c) with parameters (e, £, £)=(0.16, 
10, 0.20) as the reference model in this paper. 



ties in £ introduced by the uncertainties in the intrinsic 
scatter is not significant if this uncertainty in the scatter 
is less than 0.1 dex (for example, it is about 0.06 dex in 
Tundo et al. 2007). 

In the above models, we adopt the local BHMF es- 
timated from the M. o ~ & relation and the velocity- 
dispersion distribution function. Arguably the Af #i o — 
ibuige relation may be favored, at least for the most mas- 
sive galaxies (see Lauer et al. 2007a, but Batcheldor et al. 
2007 and Graham 2008). 5 If we adopt the local BHMF 
estimated from the M, y0 — ibuige relation and the galaxy 
luminosity function without correction of the bias as dis- 
cussed in § 2, the time integral of XAGN LF can be 
matched by Tj?°^ AGN inferred from the local BHMF if 

£ = 10, C ~ 0.15 - 0.3 but e ~ 0.08 ± 0.02±° 03 (corre- 
spondingly the spin parameter a is in the range from 0.1 
to 0.8); and therefore the QSO lifetime constrained here 
is smaller than that constrained by the local BHMF ob- 
tained from the M #i o — <J relation by a factor ~ 2. This 
is primarily due to the fact that the shape of the local 

5 Although the M. o — ^bulge relation may be favored according 
to the observations for BCGs, which primarily infer more massive 
BHs in BCGs compared with that inferred from the Af.,o — v re- 
lation, the most massive BHs may mostly be found in galaxies less 
massive than BCGs if the intrinsic scatter in the M. : o — a (or 
M. o — £<bulge) relation is significant (e.g., > 0.3 dex). 



BHMF estimated from the M,^q — £buigc relation is sim- 
ilar to that estimated from the M, t o — <J relation except 
that the normalization differs by a factor close to 2 (see 
discussions in § 2 and the bottom left panel of Fig. 1). 

5. CLUES ON THE LUMINOSITY DEPENDENCE OF THE 
OBSCURATION FRACTION OF AGNS 

Many X-ray studies have shown that the fraction of 
type 2 (or heavily obscured) AGNs decreases with in- 
creasing X-ray luminosity (Ueda et al. 2003; Akylas et 
al. 2006; Muller & Hasinger 2007, as shown by the red 
solid line and data points in Fig. 8), though there are still 
some uncertainties about whether this relation is real or 
just a selection effect (e.g., La Franca et al. 2005; Treister 
& Urry 2006; Akylas et al. 2006; Tozzi et al. 2006). This 
relation may be explained in the current evolutionary 
model for QSOs/AGNs (e.g., Hopkins et al. 2005), i.e., 
QSOs/AGNs in their rapid growth phase are moderately 
luminous and more likely to be heavily obscured, and 
as the AGN luminosity increases the UV- X-ray photons 
emitted from the QSOs/AGNs destroy the surrounding 
absorbing material and the QSOs/AGNs become unob- 
scured. 

In § 4, we have shown that the time integrals of 
the X-ray LF estimated from observations can be well 
matched by those inferred from the local BHMF within 
the reference model of the growth of individual MBHs. 
In the reference model [i.e., model (c) with parameters 
(e,f,C) = (0.16,10,0.20)], however, type 1 and type 2 
AGNs are not distinguished. In this section, we use a 
simple toy model to check whether the dependence of 
the fraction of type 2 AGNs on the X-ray luminosity can 
be really due to evolutionary effect described in the pre- 
ceding paragraph. 

In our toy model, we assume that those QSOs/AGNs 
at their early rapid growth stage are all obscured while 
those in their late evolutionary stage are all un-obscured. 
The inferred fraction of obscured QSOs/AGNs is shown 
in Figure 8 if all QSOs/AGNs in the first rapid growth 
phase (i.e., the Eddington accretion stage) with r in 
the ranges [0, r P - r s ], [0, r P - 0.5r s ], [0, r P - 0.1r s ], 
or [0, Tp] are assumed to be obscured (shown as long- 
dashed, short-dashed, dotted and dot-dashed lines in 
Fig. 8). The observations of the fraction of type 2 AGNs 
as a function of the X-ray luminosity is also shown in 
Figure 8. Although the short-dashed line in Figure 8 is 
not inconsistent with the observational trend at the high 
luminosity end [Lx ^ 10 erg s _1 ), clearly the trend 
of the dependence of the fraction of type 2 AGNs on the 
X-ray luminosity implied by these toy models is in con- 
tradiction with observations at the low-luminosity end, 
which suggests that the obscuration of AGNs cannot be 
solely an evolutionary effect arising from their individ- 
ual evolution after their nuclear activities are triggered. 
Some other effect, such as those introduced by the reced- 
ing torus model (e.g., Lawrence 1991; Simpson 2005) in 
which the opening angle of the torus is smaller in less lu- 
minous QSOs/AGNs, should be responsible for the larger 
fraction of type 2 AGNs at low luminosity. 

Comparison of the time-integral of the QSO LF in the 
optical band inferred from the local BHMF with that 
from observation will also provide information on the 
fraction of obscured QSOs in the optical band. In the 
upper panel of Figure 9, the black line represents the in- 
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ferred time integrals of the QSO LF in the B band for 
the reference model. Correspondingly, the inferred time 
integral in units of the time integral obtained from the 
QSO LF given by Richards et al. (2005) is shown in the 
middle panel. The reference model well matches the ob- 
servations at Mb < —28 but predicts more QSOs than 
those observed at magnitude Mb > — 28, which suggests 
that there exist a larger fraction of optically obscured 
QSOs/AGNs and this fraction is shown in the bottom 
panel. The dependence of the optically obscured QSO 
fraction on the luminosity in the range —27 < Mb < —20 
is much weaker compared to that in the X-ray band as 
shown in Figure 8. As we can see from the bottom panel, 
the fraction of optically obscured QSOs/AGNs can be 
as high as 80% at M B 20 23 and slightly de- 
creases to 60% at M B 27. The fraction of - 80% 

at Mb ~ —20 23 is consistent with the observa- 
tions that the ratio of Seyfert 2 galaxies to Seyfert 1 
galaxies is about 4:1 in the nearby universe. The frac- 
tion of - 60 - 70% at M B 24 27 is consistent 

with the latest estimates from Reyes et al. (2008) as in- 
dicated by the two lower limits, which are converted from 
the fractions at the [OIII] 5008A luminosity measured in 
Reyes et al. (2008) to the fractions at the B-band mag- 
nitude, and this consistence supports the constraints on 
the growth of MBHs obtained above by applying the ex- 
tended Soltan argument to the X-ray data. At higher 
luminosities, Mb < —28, the fraction of optically ob- 
scured QSOs sharply decreases to 0, which may be not 
genuine but due to effects of uncertainties in the BC at 
the high-luminosity end or the local BHMF at the high- 
mass end. 

6. THE EDDINGTON-RATIO DISTRIBUTION IN QSOS 

Observational determination of the Eddington ratio 
distribution in QSOs can put additional constraints on 
MBH growth. These are independent of, but should 
be consistent with, the constraints obtained above from 
the extended Soltan argument. Recent observational ad- 
vances allow us to seriously estimate the Eddington ra- 
tio distribution in large samples of QSOs. For example, 
using the virial mass estimators Kollmeier et al. (2006) 
and Shen et al. (2007) have shown that the logarithm of 
Eddington-ratio distribution in high-luminosity QSOs re- 
sembles a Gaussian distribution with mean around 10 -06 
to 10 -11 and width typically of 0.3 dex (see also Netzer 
et al. 2007), which may suggest that MBHs obtain most 
of their mass through accretion with a rate close to the 
Eddington limit. In this section, we check whether the lu- 
minosity evolution of individual QSOs constrained above 
is consistent with the observational Eddington ratio dis- 
tribution. 

In § 4, we have shown that the time-integrals of XAGN 
LF inferred from the local BHMF can be well matched 
to the observations using model (c) with parameters (e, 
£, C) = ((X16, 10, 0.20) for the luminosity evolution of in- 
dividual QSOs. In this reference model, the luminos- 
ity evolution and correspondingly the Eddington ratio 
evolution of a QSO arc illustrated in Figure 10. As 
shown in the upper panel of Figure 10, the luminosity 
of a QSO exponentially increases to its peak luminos- 
ity with the Eddington rate set by the self-regulation of 
disk accretion when the fuel is over-supplied, and then 
decays with time as a power-law set by the self-similar 
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Fig. 8. — Fraction of type 2 AGNs as a function of X-ray luminos- 
ity. Circles with error bars are the observed fraction of obscured 
AGNs and the solid line is the best-fit model to the red circles 
given by Akylas et al. (2006). Other lines represent the fraction 
of obscured AGNs expected from the reference model if we assume 
that all AGNs in the first rapid growth phase with r in the ranges 
[0, Tp — Tg] (long-dashed line), [0, rp — 0.5rg] (short-dashed line), 
[0, rp — O.Its] (dotted line), or [0,rp] (dot-dashed line), are ob- 
scured, while other AGNs including those in the second (declining) 
phase are all unobscured. The dot-dashed line seems to be consis- 
tent with observations at high luminosities (Lx ^ 10 43,5 erg s _1 ), 
but all blue lines arc significantly lower than the observations at the 
low luminosity end, which suggests that the obscuration of AGNs 
cannot be solely an evolutionary effect arising from their individual 
evolution after nuclear activity is triggered. 

evolution of disk accretion when the fuel is substantially 
under-supplied. The period for the QSO to have lumi- 
nosity larger than 10% of its peak luminosity is only a 
few times the Salpeter timescale. Correspondingly the 
Eddington ratio of the QSO is initially about 1 and then 
also decays with time approximately as a power-law (the 
bottom panel of Fig. 10). The timescale for m declining 
from 1 to 0.1 is relatively short compared to the Salpeter 
timescale, and those QSOs around its peak luminosity 
should mainly accrete material via Eddington ratio close 
to 1. As mentioned in § 1, the QSO LF at different red- 
shifts involves the dependence on both the nuclear ac- 
tivity triggering rate (J(z;M..o) and the luminosity evo- 
lution of individual QSOs £(t; M. i0 ) after their nuclear 
activity being triggered, so does the QSO Eddington- 
ratio distribution at different rcdshifts. Below we de- 
fine a 'time-integrated' Eddington-ratio distribution in 
QSOs, which only involves the accretion-rate evolution 
of individual QSOs. 

With a luminosity evolution model, the true Eddington 
ratio (m r ) distribution at a fixed bolometric luminosity 
Lbo\ for a MBH with present-day or final mass M. i0 is 

P(m r |L bo i, M., )dm r = (5(m r - m )dm T , (46) 

where mo = 1 if the luminosity evolution is in the first 
rapidly increasing phase, m = M^Edd^boiV-M^ if in 
the decline phase, and Ml is the mass of the MBH in 
a QSO with bolometric luminosity Lboi and with its fi- 
nal mass M. o- The M^~ can be directly obtained from 
equation (40) or (41) for fixed Lboi and M,, - With given 
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Fig. 9. — Model (c) for the B-band. The upper panel repre- 
sents the comparison of the time integral of the B-band QSO LF 
and that inferred from the local BHMF by adopting the reference 
model. The black line represents the value inferred from the local 
BHMF by using the reference model, i.e., model (c) with param- 
eters (e,£,£) = (0.16,10,0.20) (similarly in middle and bottom 
panels). The magenta line represents the values obtained from the 
fitting-formula of the QSO LF obtained by Richards et al. (2005). 
The points have similar meanings as in Fig. 3. The middle panel 
shows the ratio of the values represented by the black line and the 
points in the upper panel to the values represented by the ma- 
genta line, i.e., the values in units of that inferred from the fitting- 
formula of the QSO LF. The bottom panel shows the fraction of 
optically obscured AGNs inferred from the reference model. The 
upper arrows show the lower limit of the fraction of optically ob- 
scured AGN at z < 0.3 (M B ~ -22.9 26.1) and 0.3 < z < 0.83 

(M B ~ -26.1 27.5), respectively (Reyes et al. 2008), which are 

consistent with the prediction of the reference model in this paper. 
The range of is converted from the luminosity range of [OIII] 
5008A line [using log(L[OIII]/L e ) = -0.38A/ 2 400 -0.62 in Fig. 11 
in Reyes et al. (2008), where M2400 is the absolute magnitude at 
2400A in the rest frame, and Mb ~ Af2400 + 0.13 by assuming a 
canonical optical spectral slope (~ 0.5) of QSOs]. 



P(m r |Lboi, M,fi), the ('time-integrated') probability dis- 
tribution of the Eddington ratios among QSOs at a given 
Lboi can be defined by 



P(m r |L bo i) = / n(M. )0 )Ti t (M.,o) x 
P{L hol \M. fi )P(m r \L hoh M., )dM. )0 x 

1 -1 

n(M.,o)Tit(M.,o)P(£boi|M.,o)dM.,o 




(47) 



Note that the denominator in the above equation is just 
the time integral of the QSO LF. 



Fig. 10. — The luminosity and Eddington ratio (i.e., the accre- 
tion rate in units of the Eddington rate) evolution curves in the 
reference model, i.e., model (c) with parameters (e, £, C) = (0-16, 10, 
0.20). In the initial rapid accretion phase, the luminosity of a QSO 
exponentially increases and the Eddington ratio is a constant (~ 1) 
due to the self-regulation of the disk accretion when the accretion 
material is over-supplied, and both the luminosity and Eddington 
ratio of the QSO are followed by a rapid power-law-like decline 
due to the exhaustion of fuel and the self-similar evolution of disk 
accretion (see also Yu et al. 2005). 



Adopting the reference model for the luminosity evo- 
lution of individual QSOs, we calculate the probability 
distribution of underlying Eddington ratios among QSOs 
at a given Lboi- As shown in Figure 11, the probabil- 
ity of finding objects with low Eddington ratios in low- 
luminosity QSOs is larger than that in high-luminosity 
QSOs. The average Eddington ratio in high-luminosity 
QSOs is larger than that in low-luminosity QSOs and 
the width of the Eddington-ratio distribution in high- 
luminosity QSOs is narrower than that in low luminosity 
QSOs, although the Eddington-ratio (or the accretion- 
rate) evolution in individual QSOs is assumed to be uni- 
form. The 5-function like distribution at m r = 1 [where 
we use -^j^ exp(— x 2 /a?) with a — 0.1 to mimic the Dirac 

function 5(rh T — 1) for convenience] for each given Lboi 
represents the self-regulated rapid accretion phase with 
a rate close to the Eddington limit when the accretion 
material is over-supplied. Given an Lboi, a lower Ed- 
dington ratio corresponds to a higher MBH mass, and 
the exponential decline of the probability distribution 
of Eddington ratios at small rnJ for QSOs at a given 
Lbo\ is primarily due to the exponential-like decay of 
MBH abundance at the high-mass end (M, > 10 8 M Q ). 
The Eddington ratios in most of the luminous QSOs 
{Lho\ > 10 45 7 erg s _1 ) are close to 1 because of the 
steep falloff of the BHMF at the high-mass end and the 
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rapid decay of Eddington ratios with time in the declin- 
ing phase of the accretion-rate evolution in individual 
QSOs. These underlying Eddington-ratio distributions 
are clearly different from those observational estimates, 
i.e., a Gaussian distribution of m obs with peaks around 
10-0-6 _ iQ-i i (Kollmeier et al. 2006; Shen et al. 2007; 
Netzer et al. 2007). 

The observationally estimated Eddington-ratio distri- 
bution may be biased from the underlying true Edding- 
ton ratio distribution. The reasons are: (l) the masses of 
MBHs in QSOs are usually obtained by using the virial 
mass estimator (s) M™, and the virial mass estimator 
is based on the analysis of broad emission line reverber- 
ation mapping data for several tens of low-luminosity 
AGNs at low redshift and a calibration of it to the local 
M. o — cr relation. The estimates of M™ may scatter 
around and be offset from the real Ml, as the relation 
between luminosity and broad line region (BLR) size and 
the relation between FWHM of emission lines and BLR 
virial velocity, adopted in the virial mass estimator(s), 
are not perfect, and its validation for high luminosity 
QSOs at high redshift is not fully tested (e.g., Kaspi et 
al. 2007). A scatter of 0.3 dex in inferred M™ is plau- 
sible as pointed out by Kollmeier et al. (2006) (see also 
Shen et al. 2007) because the relation between observed 
line width and MBH mass may depend on the viewing 
angle of BLR (e.g., Krolik 2001) and the relation between 
BLR size and luminosity has an intrinsic scatter about 
0.1 - 0.2 dex (Kaspi et al. 2005). (2) There may be some 
systematic errors as large as a factor of 3 or more either 
up or down in the virial mass estimator(s) due to various 
effects, such as, a broad radial emissivity distribution, 
and an unknown angular radiation pattern of line emis- 
sion (see Krolik 2001). These systematic errors may in- 
troduce an offset of the virial mass estimator(s) from the 
underlying true mass. (3) The bolomctric luminosities 
are usually obtained using a uniform bolometric correc- 
tion (see Kollmeier et al. 2006; Shen et al. 2007), but the 
real bolometric corrections may scatter around this uni- 
form mean value by ~ 0.1 dex in the optical band (see 
§ 3.3). The dominant bias is probably those introduced 
by the viral mass estimator (s). 

We assume that the probability distribution of M™ 
for a given underlying real MBH mass Ml is 



exp 



P(logM. vir |logM. r ) 



1 



exp 



v27rAi ogM vir 

(log Af™ - logM; - e logM H 2 
2A 2 



(48) 



where A log M vi r is the scatter of MBH masses estimated 
by using the virial mass estimator(s) around the un- 
derlying given true mass, and logM vi r is the offset of 

logAf™ from the true mass of MBHs log Ml. As dis- 
cussed above, it is plausible that Ai ogM vir ~ 0.3 dex and 

|@i og Mvir| < (0.3 - 0.6) (e.g., Krolik 2001; Kollmeier et 
al. 2006; Shen et al. 2007). For a given m r at fixed iboi 
and M >i0 , the observationally estimated Eddington ratio 
is thus given by 



P(m ohs \m T ) 



(log m obs — log m r + ©i, 
2A 5 



(49) 



Combining this probability distribution with equation 
(47), the 'time- integrated' observational Eddington ratio 
(rh obs ) distribution can be inferred from the local BHMF 

as 



P(m obs |L bol ) = J dm r P(m obs \m r ) J n(M., ) X 
n t (M.,o)P(Lboi|M.,o)P(m r |£boi, M.,o)dM., x 



/ 



n(M. i o)7it(M. i o)P(Lboi|M. i o)dM. i0 



(50) 



27rln(10)A log M v.r m 



obs 



provided that the accretion rate or luminosity evo- 
lution of individual QSOs, i.e., £(t;M.. ) and thus 
P(-Lb o i|M #;0 ), is known. 

We show P(m obs |L bo i) calculated from the reference 
model in Figure 12, using equation (50) and assum- 
ing A logM vir = 0.3 dex and 6 logM vir = 0-6 dex (for 
which the Eddington ratio in the first rapid accre- 
tion phase is m r = m p = 1). It appears that the 
'time-integrated' Eddington-ratio distribution is approx- 
imately a Gaussian distribution at any fixed bolomctric 
luminosity i bo i > 10 45 75 erg s _1 but with a small tail 
at the low-Eddington ratio end. The Gaussian-like dis- 
tribution mainly corresponds to the peaks at m r = 1 
as shown in Figure 11, which represent the rapid accre- 
tion phase of individual QSOs with a rate self-regulated 
by the Eddington limit. The width of the Gaussian- 
like distribution mainly reflects the scatter A logM vir in 
the estimates of MBH masses using the virial mass es- 
timator^), and the locations of peaks in the Eddington 
ratio distribution arc roughly determined by the offset 
©log M vir and the value of Eddington ratio m P during the 
self-regulated rapid accretion phase when the accretion 
material is over-supplied. For QSOs with lower bolomct- 
ric luminosities (Lboi S 10 45 25 erg s _1 ), the probabil- 
ity of finding low Eddington-ratio (m obs < 0.03) objects 
becomes significant, which is primarily because the un- 
derlying MBH mass function is shallow at the low-mass 
end (M. < 10 8 M Q ) and the decline phase of the self- 
similar evolution of the disk accretion (see also Fig. 10) 
around big MBHs contributes significantly to the counts 
of low bolometric luminosity objects. Therefore, the Ed- 
dington ratio distribution among low-luminosity QSOs 
should provide independent constraints on the long-term 
evolution of disk accretion, especially in the decline phase 
(see also Yu et al. 2005). 

Except for giving a rough comparison with observa- 
tions below, we do not intend to use P(m obs |L bo i) in- 
ferred from the reference model in this paper to directly 
fit the observational Eddington ratio distribution esti- 
mated by Kollmeier et al. (2006), Shen et al. (2007), and 
Netzer et al. (2007). The reasons are: (1) P(m obs |L b oi) 
obtained from equation (50) is a 'time-integrated' and 
volume-weighted Eddington-ratio distribution, while the 
current observationally estimated distributions are for 
given redshift intervals and not volume weighted. If 
the Eddington-ratio distribution at a given bolometric 
luminosity is independent of redshift as suggested by 
Kollmeier et al. (2006) (but perhaps it is not as argued 
by Shen et al. 2007), the 'time- integrated' Eddington ra- 
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Fig. 11. — The underlying 'time-integrated' Eddington ratio dis- 
tribution among QSOs at a given bolomctric luminosity inferred 
from the reference model in § 4. From left to right at the low- 
Eddington ratio end, the lines represent P(m r \L\ 10 i) at the given 
bolometric luminosity -Lboi = 10 43 ' 75 (solid line), 10 44 ' 25 (dot- 
ted line), 10 44 ' 75 (short-dashed line), 10 45 - 25 (long-dashed line), 
10 45.75 (dot-short-dashed line), 10 46 ' 25 (dot-long-dashed line), 
iq46.75 ( s hort-dash-long-dashed line), and 10 47 ' 25 erg s _1 (solid 
line), respectively. See details in § 6. 

tio distribution among QSOs may represent the obser- 
vational one for different redshift intervals; (2) the ob- 
servational Eddington ratio distribution may be biased 
significantly at low Eddington ratios in flux-limited sur- 
veys, such as the SDSS (sec Kollmeier et al. 2006 and 
Shcn ct al. 2007; for more general discussion of selection 
bias see also Lauer et al. 2007b and Yu & Lu 2004b); (3) 
P(m obs |Lboi) obtained from equation (50) does not dis- 
tinguish obscured and unobscured QSOs, while the ob- 
servationally estimated distribution is primarily obtained 
from optical QSO samples. If the obscuration is only a 
geometrical effect, then P(m obs |Lboi) may be the same 
as the Eddington-ratio distribution in optical QSO sam- 
ples; however, if the obscuration is partly due to an evo- 
lutionary effect (e.g., QSOs may be more likely obscured 
in their early accretion stage), the probability of QSOs 
with high Eddington ratios may be suppressed. 

If the MBH masses inferred from the virial mass esti- 
mator have an offset ~ 0.5 — 0.6 dex from the underly- 
ing true masses and also a scatter of ~ 0.3 dex around 
them, P(m obs |Lboi) at any fixed bolometric luminosi- 
ties Lboi ^ 10 45 75 erg s _1 is roughly consistent with the 
Gaussian-like distribution of Eddington ratios (with a 
scatter of ~ 0.3 dex, and peak at 1/4) among QSOs esti- 
mated by Kollmeier et al. (2006); however, P(m obs |Lboi) 
at bolometric luminosity Lboi ^ 10 45 25 erg s" 1 has a tail 
extended to lower Eddington ratios (m obs < 0.01), which 
appears to conflict with that obtained by Kollmeier ct al. 
(2006) for QSOs with L bol < 10 45 5 erg s" 1 . The Edding- 
ton ratio in the self-regulated rapid accretion phase is 
assumed to be 1 (i.e., exactly the Eddington limit) in the 
above calculations, but it may be slightly different from 
1 and a little smaller. For example, if rh 1 — rh p = 0.5 
in the first rapid accretion phase of the model (c), the 



time- integral of XAGN LF inferred from the local BHMF 
can still marginally match that obtained from observa- 
tions 6 , but it cannot if m r = m p is substantially less 
than 0.5. If m P — 0.5, to reproduce the observational 
Eddington ratio distribution at high bolometric luminos- 
ity (L bo i > 10 45 75 erg s" 1 ), O - 0.2-0.3 dex is required; 
and the tail of P(m obs |Lboi) at low Eddington ratios for 
low-luminosity QSOs (L bo i < 10 45 25 erg s _1 ) is still sig- 
nificant, though less significant compared to the case for 
m r = Trip = 1 in the first rapid accretion phase. 

The consistency between P(m obs |Lboi) inferred from 
the extended Soltan argument and the observationally 
estimated Eddington-ratio distribution at high bolomet- 
ric luminosity (Lboi ^ 10 45 75 erg s~ x ) suggests that the 
majority of bright QSOs accrete material at a single rate 
close to the Eddington limit. But we should be cautious 
of any over-interpretation of the possible inconsistency at 
low bolometric luminosity (Lboi ^ 10 45 25 erg s _1 ) above, 
since the observational results obtained by different au- 
thors have not yet converged (e.g., Kollmeier et al. 2006; 
Shcn ct al. 2007; Netzer et al. 2007). For example, the 
observational Eddington-ratio distribution among QSOs 
with L bo] < 10 45 5 erg s" 1 obtained by Shen et al. (2007) 
(i.e., a Gaussian-like distribution with peak at 10 -1 ' 1 and 
a scatter of 0.42 dex) is substantially shifted to lower Ed- 
dington ratios compared with that obtained by Kollmeier 
et al. (2006) (i.e., a Gaussian-like distribution with peak 
10~ - 6 and scatter 0.3 dex). Note also that the Eddington 
ratios in low-luminosity AGNs at low redshifts do cover a 
wide range as shown by Woo & Urry (2002), Heckman et 
al. (2004), and Greene & Ho (2007), which may be con- 
sistent with the prediction obtained from equation (50) 
above. 

The offset logM vir <~ 0.3 — 0.6 dex required by the 
above observational constraints suggests that the MBH 
masses inferred from the virial mass estimator(s) may 
be over-estimated by a factor of 2 — 4, at least for high- 
luminosity QSOs, which is compatible with the possi- 
ble systematic errors in the MBH mass estimated by the 
reverberation mapping technique (Krolik 2001). If this 
offset is real, it is intriguing and important since many 
current studies on the growth and evolution of MBHs 
in QSOs are based on the virial mass estimator(s). For 
example, the masses of MBHs in two samples of AGNs 
at redshifts z = 0.36 and z = 0.57, estimated from the 
virial mass estimator(s), are found to be larger than that 
estimated from the local M #! o — o relation by 0.54 dex 
and 0.51 dex, respectively (Treu et al. 2004; Woo et al. 
2006, 2008), which is suggested as an indicator of that 
the growth of MBHs predates the final growth of bulges 
in these AGN host galaxies. If MBH masses from the 
virial mass estimator(s) are generally over-estimated by 
a factor of 2 — 4 as argued above, then there should be 
not much difference between the rescaled virial masses 
and that predicted from the M, — a relation for those 
MBHs in the studies of Treu et al. (2004) and Woo et al. 
(2006, 2008). 

6 Our calculations show that this set of the accretion rate in the 
first rapid accretion phase seems to under-predict the time-integral 
of the optical QSO LF at the high-luminosity end (Mb < —28), 
which may be partly due to some uncertainties in the bolometric 
correction for the B band at the high-luminosity end; otherwise 
QSOs must accrete at a rate closer to the Eddington limit during 
the first rapid accretion phase. 
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From left to right at the left side of the lines: 
L bol = 10 43 - 75 to 10 47S5 erg/s, AlogL M =0.5dex 
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Fig. 12. — The observational 'time-integrated' Eddington ratio 
distribution among QSOs at a given bolometric luminosity inferred 
from the reference model in § 4 by assuming that the scatter and 
offset in the masses of MBHs obtained from the virial mass esti- 
mator^) are 0.3 dex and 0.6 dex, respectively. The line types have 
the same meaning as in Fig. 11. 



Similar to the probability distribution of the Edding- 
ton ratio among QSOs at a given bolometric luminosity 
given in equation (47) , the probability distribution of the 
Eddington ratio among QSOs at a given MBH mass Ml 
can be estimated by 

P(m obs |M. r ) = J dm r P(m obs |m r ) / n(M. i0 ) x 
nt{M. fl )P{Ml\M. fl )P{m r \Ml M. fi )dM. fl x 



n(M. fi )n t {M. fi )P(M r .\M. fi )dM. fi 



(51) 



or the probability distribution of Eddington ratio among 
QSOs at a given M™ can be estimated as 

P(m obs |M. vir ) = J drn?P{m ohs \rrf) J dM. r P(M. vir |M. r ) x 
n{M. fi ) nt {M. fl )P{Ml\M. fl )P{m r \Ml,M. fl )dM. fl / 

J P{M™\Ml)dMl J n(M.,o)r lt (M. ; o)P(M. r |M. ; o)dM. ;0 , 

(52) 

where P(M. vir |M. r ) = P(logM™| log M. r )/[M. vir ln(10)]. 
We note here that P(rn obs |M^ ir ) should be skewed to- 
ward low Eddington ratios because a QSO may spend a 
majority of its lifetime in the declining phase with small 
m r (< 0.1) as that revealed by the reference model in this 
paper (e.g., Fig 10). However, it is not easy to obser- 
vationally estimate P(m ohs \M™) since low-luminosity 
QSOs are more likely to be missed in flux-limited sur- 
veys, especially at high redshift (see also discussions in 
Kollmeier et al. 2006). Wc defer the comparison of this 
distribution with observations to future work. 



7. TOY MODELS FOR MERGERS 
7.1. Gas-poor (dry) mergers 

In current hierarchical galaxy formation models, merg- 
ers of galaxies are the main route to form elliptical galax- 
ies and stellar bulges (e.g., Kauffmann, White & Guider- 
doni 1993; Cole et al. 2000). If each merging galaxy has a 
central MBH, mergers of two galaxies will inevitably form 
binary MBHs and may further lead to mergers of MBHs if 
their inspiral and orbital evolution time is shorter than a 
Hubble time (e.g., Begelman, Blandford & Rees 1980; Yu 
2002). Mergers of MBHs occurred after the quenching of 
nuclear activity, as the probable consequence of galaxy 
gas-poor (dry) mergers, may significantly re-shape the 
BHMF established by accretion processes. Below we use 
a toy model to illustrate the change of the BHMF due to 
BH mergers under the assumption that MBHs grow only 
by merging BHs but with little accretion directly onto 
MBHs during the galaxy dry merger stage. 

Assume that a major dry merger of two host galax- 
ies always leads to the merger of their central MBHs 
with mass ratio of a = M. ; 2/(M. i i + M. j2 ) (e.g., 0.5 or 
0.25 for a 1:1 or 1:3 merger) after the nuclear activity 
is quenched, where M #j i and M #j 2 are the masses of the 
two BHs before their merger. The merged BH mass is 
M. i0 = /3(M.,i + M. ;2 ) and (1 - (3) is the fraction of en- 
ergy (or mass) losses due to gravitational waves during 
the MBH merger process. Thus, the mass function of 
QSO remnants right after the nuclear activity is 



n' M .( M ',,o) = J rfM. ;0 n M .(M. ; o) x 



S(M'. fl - -M., ) + S(M' mfi 



: M., ) 



= -n M .( M ;o)\ M .,, 



+ 



-n M .(M., )\ M j!_ 



1 _ a ..,„. v -.,„ / i M , =J-M: i0 ' ( 53 ) 

Note this mass function n' M% {M', a ) is non-synchronous 
since the last major (dry) merger may occur at different 
time for different MBHs. 

During the merging process of BHs (which is di- 
vided into three phases: inspiral, merger and ringdown), 
the total energy lost through gravitational waves, E ra d 
is difficult to calculate, especially for the merger of 
two BHs with large spins, but roughly in the range 
0.03M., 12 F( M /M M2 ) < £ rad < 0.2M., 12 F( M /M., 12 ), 
where M #l i2 = M #i i + M #i 2 is the total (initial) mass of 
the two BHs, \i is the reduced mass and F(/x/M #i i 2 ) = 
(4fi/M.. 12 ) 2 (see eq. 3.7 in Flanagan & Hughes 1998). 
Recent breakthrough in relativistic numerical calculation 
of merging binary BHs due to Pretorius (2005) and Baker 
et al. (2006) has shown that on the order of 5% of the 
initial rest mass for a system of two equal mass, nonspin- 
ning BHs is radiated as gravitational waves during the 
final orbit and ringdown, which is consistent with the es- 
timate in Flanagan & Hughes (1998). For the merging of 
equal mass, rapidly spinning BHs, the energy radiated as 
gravitational wave could be larger and the upper limit is 
about 24% if the final spin is around 0.9 (e.g., Pretorius 
& Khurana 2007). Since the spin of MBHs is probably 
close to 0.7 — 0.9 due to accretion processes (e.g., Gam- 
mic et al. 2004; Shapiro 2005; Hawley et al. 2007), here 
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we choose two cases, 10% and 24% of the initial total 
rest mass, for the amount of energy radiated as gravita- 
tional waves, and therefore (3=1 — 0.1F(/z/M #j i 2 ) and 
1 - 0.24F(/i/M. i i 2 ), respectively. 

Using the data from Galaxy Evolution from Morphol- 
ogy and SEDs (GEMS), Bell et al. (2006) find that 
present-day spheroidal galaxies with My < —20.5 on 
average have undergone between 0.5 and 2 major dry 
mergers since redshift z < 0.7 (see also similar results 
in Conselice et al. 2003 for redshift z < 3, and Lin et al. 
2004 for redshift z < 1.2). MBHs with mass substantially 
less than 10 8 M© are mostly hosted by the stellar bulges 
of spiral galaxies, which should not have undergone a sig- 
nificant number of major dry mergers in the near past 
(e.g., z < 1) since their disks are preserved. For simplic- 
ity, here we assume all MBHs with M. j0 > 1O 8 M , cor- 
responding to My < —20.5, have experienced one major 
dry merger after the quenching of nuclear activity, while 
smaller MBHs did not experience major dry mergers. Us- 
ing equation (53) to correct the effect due to dry mergers 
in riM. (M..o), the BHMF right after the quenching of nu- 
clear activities, nM.(M. ), which is established by the 
accretion process, is estimated and shown in Figure 13. 
As seen from Figure 13, the abundance of MBHs with 
mass larger than a few 10 9 M Q may be enhanced at most 
by a factor ^2 — 3 due to the major dry mergers after 
the quenching of nuclear activities. This enhancement is 
not so prominent compared with that in the estimate of 
BHMF due to the uncertainty in the intrinsic scatters in 
the M,, — v an d M #i n — -^buige relationships. For ex- 
ample, a slight error in the estimated intrinsic scatters, 
e.g., 0.05 dex, could introduce an uncertainty larger than 
a factor of ^ 2 — 3 at the high-mass end of the BHMF. 
Therefore, we conclude that it is safe to neglect the un- 
certainty in the BHMF due to major dry mergers after 
the quenching of nuclear activities in our calculations as 
the intrinsic scatters in the M #i o - c or M,,o — £t>uige 
relations are currently not well determined. 

7.2. Gas-rich (wet) mergers 

For a present-day MBH with mass M #: o, its host 
galaxy may have undergone more than one (wet) major 
merger. By 'wet' major merger here we mean gas-rich 
major mergers (spiral+spiral or spiral+clliptical) which 
lead to substantial gas fueling to trigger nuclear activ- 
ity. In the current scenario of hierarchical formation and 
co-evolution of galaxies and MBHs (e.g., Kauffmann & 
Haehnelt 2000; Bower et al. 2006; Croton et al. 2006; 
Malbon et al. 2007), each 'wet' major merger leads to 
a nuclear active phase with substantial increase in the 
central MBH mass. For a local MBH we see today, it 
may build up through multiple times of nuclear activ- 
ity and mergers of smaller MBHs, and thus our assump- 
tion of a single time of nuclear activity in § 4 may be 
an over-simplification. In order to estimate the effects 
of possible multiple phases of nuclear activity in the as- 
sembly history of a MBH, we assume that all MBHs ex- 
perienced two or more 'wet' major mergers and its mass 
increases by a factor of 5, 10 or more after each 'wet' ma- 
jor merger, and this assumption is compatible with the 
current co-evolution models (e.g. Kauffmann & Haehnelt 
2000; Bower et al. 2006; Croton et al. 2006; Malbon et al. 
2007). Before each 'wet' major merger, we assume that 
the MBH triggered by the 'wet' merger later is the merger 
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Fig. 13. — Effect of dry mergers on re-shaping the BHMF 
after the quenching of nuclear activity. The solid line represent 
the BHMF M.^M. (M»,o) in the local universe estimated by the 
M,,o — o relation and the velocity-dispersion distribution as that 
shown in Fig. 1. The dashed and dotted lines show the BHMF 
bm, (M', ) after correcting the effect due to dry mergers. Dotted (or 
dashed) lines represent the case that all MBHs with mass > 10 s Mq 
experienced one 1:1 (or 1:3) dry major merger, while other MBHs 
with smaller mass have not undergone any major mergers, after 
the quenching of nuclear activity. The loss of energy through grav- 
itational waves are assumed to be 10% (thick dotted or dashed 
lines) or 24% (thin dotted or dashed lines) of the initial total rest 
mass. This figure shows that the abundance of MBHs after the 
quenching of nuclear activity may be smaller than that estimated 
in the local universe at most by a factor of 2—3 at the high-mass 
end, but is not significantly different from the local BHMF at the 
low-mass end. This difference caused by the effect of dry mergers 
on the estimated BHMF is not so prominent compared with that 
due to the uncertainty of intrinsic scatters in the M,fi — a and 
M..o — Lbulgc relations. 



remnant of two progenitors with mass ratio either 1 : 1 
or 1 : 3. These two progenitors also experienced a similar 
period of significant mass growth before and so on and 
so forth. We can use this procedure backwards for two 
or more times to mimic multiple phases of accretion (in- 
cluding mergers). With the above assumptions, combin- 
ing equations (53) and (24), the time- integral of XAGN 
LF can be calculated. In these calculations we also adopt 
the same parameters (e, £, 7) as in the reference model 
but adjust £ in order to satisfy the assumption of a factor 
of 5, 10 or more mass increase during each 'wet' major 
merger. We find that the inferred value of the time- 
integral of XAGN LF at the low-luminosity end for the 
multiple-times nuclear activity assumption is larger than 
that obtained by assuming a single-time nuclear activity 
at most by ~ 20 — 30%, and the difference is negligible 
at the high-luminosity end. We therefore conclude that 
the assumption of a single-time nuclear activity is good 
enough, provided that the masses of MBHs substantially 
increase (say, by a factor of 5, 10 or more) during the 
nuclear activity. 
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8. FURTHER IMPLICATIONS 

8.1. BHMF at redshift z 

If the probability function P(L\M t ) is known (for ex- 
ample, as described by the reference model we obtained 
from the extended Soltan argument above), the BHMF 
at redshift z, i.e., nM.(M.fi,t z ) can be estimated by an 
equation similar to equation (1) 

n M .( M ;0,tz)nt{M. t0 ) x 



/ 

J z 



dz 



dz 



I 

Jo 

P(L\M. fi )dM. fi , (54) 
where t 7 = f°° 14^-1 dz is the cosmic time at redshift z. 

z J z I dz\ 

Equation (54) is true only if the BHMF at z is dominated 
by quiescent MBHs. 

If MBHs with final mass M..o shincd for a time 
i~\t{M m fl) with luminosity L — ALEdd(M.) and with- 
out increasing their mass significantly (M. o — M.), 
where A is a constant and LEdd(M.) is the Eddington 
luminosity (see definition in eq. 25), then P(L\M,fl) ~ 
5[L — XLEdd(M, y o)] and thus 



n M .{ M ;0,to)' 



1 



71t(M., ) 
dt 



* l(L, z)\ L=X L Edd (M. , ) X 



dz 



dz x X 



dL Edd (M.) 



dM. 



(55) 



M.=M. |( 



If these assumptions are correct, it appears that the local 
BHMF can be forced to be always consistent with the 
QSO LF by adjusting nt(M. t o) and A. Similarly, the 
BHMF at redshift z is given by 



n M .(M, fi ,t z )i 



1 



Tlt(M. ;0 

dt 



dz' 



dz' x A 



$l(L, z)\ L= x LEdd (M. i0 ) X 

dL Edd (M.) 



dM. 



(56) 



M.=M. i0 



The BHMF at redshift z, riM.(M,.o,t z ), estimated from 
equation (56) may be not accurate if the real Tit(M., ) 
is substantially longer than the Salpeter timescale since 
in this case the MBH mass was evolving rapidly dur- 
ing its active phase and the BHMF sharply decreases 
at the high-mass end, and also the assumption of P ~ 
S(L — AZvEdd) is not good. The errors in the BHMF at 
redshift z obtained from this simple approach can be esti- 
mated by comparing it with that obtained from equation 
(54). With the constraints on the luminosity evolution 
of individual QSOs obtained above, we will estimate the 
BHMF at different redshift z and check whether these 
estimates are consistent with observations in a future 
study. 

8.2. Triggering rate of nuclear activity 
The QSO LF at redshift z can be inferred as 

* L (L, z) = J rfM.,o J g( Zl - M., )5(L - C(t; M., )) x 

(57) 



dt 



dzi 



dz,. 



where t = t z — t Zi . Given ^>l(L,z) and £(t;M.. ), 
the triggering rate of nuclear activity in the mass range 
M.,o — M.,o + dM.fi, G(z; M. j0 ), can be solved from the 



above integral equation. Estimation of C?(z;M..o) is of 
fundamental importance because it is this function that 
dominates the cosmic evolution of the QSO population 
and the down-sizing nature of the formation of MBHs. 
Here the cosmic evolution of the QSO population means 
that the comoving number density of the QSO popula- 
tion brighter than a certain luminosity (or in a certain 
luminosity range) has a peak at an intermediate redshift 
(e.g., z ~ 2 — 3) and decreases at both higher and lower 
redshift (e.g., Richards et al. 2006a; Hasinger et al. 2005) 
and the down-sizing nature refers to that observationally 
the characteristic mass of MBHs in QSOs decreases with 
decreasing redshifts (e.g., Marconi et al. 2004; Merloni 
2004). Given Q{z;M. fi ) and £(t;M. j0 ), many statisti- 
cal properties of QSOs (for instance, the Eddington rate 
ratio distribution and the MBH mass distribution at dif- 
ferent redshifts in QSOs/AGNs) can be inferred. Com- 
parison of these inferred properties with those directly 
obtained from observations will further reveal details of 
the growth and evolution of MBHs and QSOs. We will 
present this in a future study. 

9. CONCLUSIONS 

In this paper, we have studied the observational con- 
straints on the growth of MBHs using the extended 
Soltan argument. In this approach, the local BHMF is 
directly connected with the time-integral of the QSO LF 
through only the luminosity evolution of individual QSOs 
(and correspondingly the accretion-rate evolution, given 
the mass-to-energy conversion efficiency), and the lumi- 
nosity evolution of individual QSOs is isolated from the 
cosmic evolution of the triggering rate of nuclear activity. 
The luminosity (or accretion-rate) evolution of individ- 
ual QSOs has an unambiguous physical definition, and 
it is different from the 'mean accretion rate' as a func- 
tion of mass and/or redshift widely used in the literature 
(e.g., Marconi et al. 2004; Shankar et al. 2004, 2007) in 
that the 'mean accretion rate' is a combined property de- 
pending on both the luminosity evolution of individual 
QSOs and the cosmic evolution of the triggering rate of 
nuclear activity. 

With recent knowledge of the relationships between 
MBH mass and host galaxy properties (cither the M. o — 
a relation or the M., — ^buigo relation) and the distri- 
bution of galaxy properties (either a or Lbuige), we es- 
timate the local BHMF. We obtain the time-integral of 
the QSO LF from recent estimates of QSO LFs in both 
optical and X-ray bands. Using the local BHMF and 
the time-integral of the QSO LF, we obtain robust con- 
straints on the luminosity (or accretion rate) evolution 
of individual QSOs and important characteristic param- 
eters describing the growth of individual MBHs, such 
as the mass-to-energy conversion efficiency and lifetime, 
which are summarized below. 

• The luminosity (or accretion rate) evolution of in- 
dividual QSOs probably involves two phases: an 
initially exponentially increasing phase set by the 
Eddington limit (i.e., C ~ -^Edd) when the infall 
material to feed the central MBHs is over-supplied; 
and then followed by a phase with power-law de- 
clining set by a self-similar long-term evolution of 
disk accretion (i.e., C oc t~ 7 and 7 ~ 1.2 — 1.3). 
With this type of luminosity evolution, the time- 
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integral of QSO LF can be well matched by that 
inferred from the local BHMF. Other simple lumi- 
nosity evolution models, such as a single Eddington 
ratio for all MBHs/QSOs or an initially exponen- 
tially increasing phase followed by an exponentially 
decay phase, cannot satisfy the extended Soltan ar- 
gument simultaneously at both the high-luminosity 
end and low-luminosity ends, and thus are ruled 
out. 

The mass-to-energy conversion efficiency e is ~ 
0.16 ± 0mt° 05 (correspondingly the spin param- 
eter a is in the range from 0.8 to 0.99) if adopting 
the local BHMF estimated from the M. o — a rela- 
tion, which is fully consistent with the theoretical 
expectations of ~ 0.10—0.20, i.e., the spin of MBHs 
in QSOs may stay at an equilibrium of <~ 0.7 — 0.9 
for most of the QSO lifetime (e.g., Gammie et al. 
2004; Shapiro 2005; Hawley et al. 2007). However, 
the efficiency e is reduced to ~ 0.08 ± 0.02l° 03 
(and correspondingly the spin parameter a is in the 
range from 0.1 to 0.8) if adopting the local BHMF 
estimated from the M. ; o — ibuigc relation, which is 
lower than but may be still marginally consistent 
with theoretical expectations. 

The lifetime of QSOs/AGNs, which depends on de- 
tailed definition of the nuclear activity or the lower 
threshold set to the active nuclear luminosity, can 
be as long as a few 10 9 yr, and the character- 
istic timescale in the luminosity increasing phase 
and transition timescale to the declining phase do 
not necessarily depend on the mass of their central 
MBHs. The period that a QSO or MBH radiating 
at a luminosity larger than 10% of its peak lumi- 
nosity is only about 2 — 3 x 10 8 yr, and during 
this period the MBH obtained most of its mass. 
If adopting the local BHMF estimated from the 
M #j o — ibuigo relation, the above values related to 
the QSO lifetime decrease by a factor <~ 2. 

For individual QSOs, the characteristic timescale 
for the luminosity (or accretion rate) to decline 
from its peak Lp to O.ILp in their second phase 
should be relatively short compared to the Salpeter 
timescale, which suggests that the material in- 
falling from a large galactic scale and deposited in 
the vicinity of MBHs be consumed by rapid accre- 
tion onto the central MBH and at the mean time 
further deposit of material can be efficiently sup- 
pressed by some mechanisms, probably the AGN 
feedback mechanism (e.g., Silk & Rees 1998; King 
2003; Murray et al. 2005; Di Matteo et al. 2005), 
on a timescale < rs- 

The majority of high-luminosity (£boi ^ 
10 45 75 erg s _1 ) QSOs accrete material via an 
almost single Eddington ratio, close to 1 and not 
smaller than half of the Eddington limit, which 
suggests that the disk accretion onto MBHs should 
be indeed self-regulated by the Eddington limit 
when the accretion material is over-supplied in the 
initial phase. 



• Low-luminosity QSOs (£boi ^ 10 45 25 erg s" 1 ) ac- 
crete material via a much wider range of Eddington 
ratios, and a significant fraction of them accrete 
material via low Eddington ratio (m < 0.1), which 
corresponds to the self-similar long-term evolution 
of disk accretion around MBHs (M, ~ r~ 7 and 
7 ~ 1.2 — 1.3) when the accretion material is under- 
supplied. 

• The Eddington ratio distribution among 
QSOs/AGNs inferred from the extended Soltan 
argument concentrates toward high Eddington 
ratios (close to 1), especially for high luminos- 
ity QSOs, which appears to conflict with that 
estimated directly from observations (with a 
mean value of <~ 10~ 06 — 10 -11 ) by using the 
virial mass estimator(s). To make these two 
distributions consistent with each other, an offset 
of 0.3 — 0.6 dex in the MBH mass estimated from 
the virial mass estimator(s) is required, at least for 
high-luminosity QSOs, which suggests that MBHs 
masses obtained from the virial mass estimator(s), 
have been systematically over-estimated by a 
factor of 2 - 4. 

• The fraction of optically obscured QSOs/AGNs in- 
ferred from the extended Soltan argument can be 

as high as 80% at M B 20 23 and slightly 

decreases to 60% at Mb = —24 27, and these 

numbers are consistent with recent observations by 
Reyes et al. (2008). The dependence of the fraction 
of type 2 AGNs on the X-ray luminosity cannot be 
solely an evolutionary effect arising from their indi- 
vidual evolution after nuclear activity is triggered 
(i.e., QSOs are more likely to be obscured in the 
early stage of the MBH growth), and some other 
effects (e.g., those introduced by the receding torus 
model; Lawrence 1991) should be responsible for 
the larger fraction of type 2 AGNs at low luminosi- 
ties (L x < 10 43 5 erg s^ 1 ). 

We estimate possible effects due to MBH mergers 
(which may re-shape the local BHMF) and multiple 
times of nuclear activity and accretion (e.g., triggered 
by multiple times of galaxy 'wet' major mergers) in the 
growth history of a MBH, and we find that these effects 
on our conclusions are insignificant, which again supports 
that the constraints obtained above are robust. 

The constraints on the luminosity evolution of indi- 
vidual QSOs obtained from the extended Soltan argu- 
ment in this paper, together with the QSO LF, can be 
further used to derive the BHMF at high redshifts and 
the cosmic evolution of the triggering rate of nuclear ac- 
tivity. These constraints and those recent estimates on 
the Eddington ratio distribution in QSOs ask for serious 
theoretical modeling of the long-term evolution of disk 
accretion around MBHs. More detailed modeling of ac- 
cretion and radiation transfer physics in the vicinity of 
MBHs may have to be involved to determine the self- 
regulation of the disk accretion (rather than the simple 
Eddington limit argument) at the initial phase with suf- 
ficient deposited accretion material. It should be one of 
the important long-term goals for theoretical studies on 
the growth of MBHs to answer questions like what de- 
termines the transition from the initial rapid accretion 
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phase with Eddington ratio close to 1 to the rapid de- 
clining phase, what shuts off the efficient accretion pro- 
cess around MBHs (probably jointly determined by an 
efficient feedback mechanism and the accretion disk vis- 
cosity), and what determines the evolution of the spin of 
MBHs. 
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